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ABSTRACT 

This  report  describes  20  investigations  in  the  fields  of  seismic 
discrimination  and  the  analysis  of  data  from  a global  seismic 
network.  These  are  grouped  as  follows:  the  estimation  of 
event  location  and  focal  depth  (4  contributions),  the  estimation 
of  the  source  moment  tensor  (5  contributions),  miscellaneous 
studies  (6  contributions),  and  computer  systems  and  data  col- 
lection (5  contributions). 
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SUMMARY 


This  is  the  twenty-ninth  Semiannual  Technical  Summary  report  describing  the  activities  of 
Lincoln  Laboratory  funded  under  Project  Vela  Uniform.  This  report  covers  the  period  1 April 
to  30  September  1978.  Project  Vela  Uniform  is  a program  of  research  into  the  discrimination 
between  earthquakes  and  nuclear  explosions  by  seismic  means.  A recent  new  emphasis  of  the 
project  is  in  the  development  of  the  data -handling  and  analysis  techniques  that  might  be  appropri- 
ate for  the  monitoring  of  a potential  Comprehensive  Test  Ban  Treaty,  presently  under  negotia- 
tion. The  objectives  of  the  Lincoln  Laboratory  program  are  to  carry  out  fundamental  research 
with  the  seismological  problems  associated  with  the  detection,  location,  and  identification  of 
earthquakes  and  nuclear  explosions,  and  to  develop  both  processing  algorithms  and  data- 
management -system  techniques  that  may  be  implemented  as  part  of  U.S.  activities  associated 
with  an  international  treaty  monitoring  network.  Much  of  this  research  involves  the  analysis 
of  digital  seismic  data  from  the  new  network  of  high-quality  seismic  stations  called  Seismic 
Research  Observatories  (SROs). 

A fundamental  problem  in  seismic -data  analysis  is  the  development  of  rapid  automated  or 
semiautomated  methods  for  the  routine  estimation  of  source  location  and  focal  depth,  and  the 
association  of  long-period  (LP)  data  with  events  determined  from  short-period  (SP)  arrivals. 

One  investigation  describes  an  attempt  to  determine  source  azimuth  and  distance  directly  from 
surface-wave  data,  for  comparison  with  SP  data.  With  further  refinement,  this  technique  may 
greatly  assist  the  association  problem.  Event  location  is  primarily  a problem  in  path  calibra- 
tion, and  two  studies  address  the  application  of  the  master-event  technique.  If  local  crustal 
structure  is  relatively  simple  and  well  known,  it  is  shown  that  the  addition  of  regional  data, 
using  travel  times  of  crustal  phases,  can  provide  a noticeable  improvement  to  the  quality  of 
epicenters  determined  from  teleseismic  data  alone.  The  master-event  technique  can  also  be 
useful  in  the  estimation  of  the  relative  focal  depth  of  events.  In  another  study,  the  application 
of  seismic  ray  tracing  to  event  location  in  a complex  area  is  described. 

We  have  begun  a series  of  investigations  into  methods  for  the  estimation  of  fault  plane  so- 
lutions and  the  moment  tensor.  Two  studies  investigate  the  direct  inversion  of  body-wave  ampli- 
tudes to  obtain  the  moment  tensor.  In  the  first,  a case  is  described  in  which  inversion  leads  to 
a predominantly  double -couple  source  which  has  an  orientation  rather  more  consistent  with  the 
tectonic  environment  than  the  published  fault  plane  solution.  The  second  study  attempts  to  es- 
timate the  potential  errors  introduced  into  the  moment  tensor  inversion  by  structural  inhomo- 
geneity. Improvements  to  the  determination  of  fault  plane  solutions  depend  critically  on  the 
development  of  nonsubjective  approaches.  One  such  approach,  using  a linear-estimation  pro- 
cedure, is  described.  Another  study  outlines  the  possible  utilization  of  synthetic  seismograms 
computed  by  the  superposition  of  normal  modes,  and  formulates  the  inverse  problem  of  esti- 
mating the  moment  tensor  from  observed  seismograms. 

A number  of  miscellaneous  studies  are  included,  consisting  mainly  of  analysis  of  SRO  data. 
Two  investigations  have  focused  on  the  problem  of  event-detection  algorithms.  The  first  as- 
sesses the  performance  of  the  current  SRO  event  detector,  and  demonstrates  a high  incidence 
of  false  alarms  and  disappointing  performance  at  some  stations.  The  second  discusses  some 
potential  improvements  to  detection  algorithms.  Another  study  has  begun  to  assess  the  useful- 
ness of  digital  waveforms  in  the  estimation  of  body-wave  magnitude  m^.  Station  amplitude  bias 


r 
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appears  to  be  a significant  effect,  and  calibration  by  source  region  may  be  possible.  A study 
of  mantle  Love-wave  dispersion  is  a beginning  of  a global  analysis  of  regional  variations  in 
phase  velocity  and  Q.  Rayleigh-wave  dispersion  is  used  to  estimate  the  structure  of  the  litho- 
sphere in  a remote  aseismic  area  of  the  South  Atlantic.  Investigations  of  the  propagation  of  the 
crustal  phase  Lg  continue.  A study  shows  that  the  observed  phase  velocities  of  the  vertical 
component  of  Lg  are  consistent  with  the  expected  behavior  of  higher-mode  Rayleigh  waves. 

We  continue  to  make  improvements  to  our  in-house  computer  system,  and  are  focusing 
our  programming  effort  into  those  areas  that  will  have  eventual  application  in  a global  treaty 
monitoring  context.  A new  version  of  our  Data  Analysis  and  Display  System  (DADS)  is  described, 
together  with  improvements  in  our  waveform  graphics  package  and  the  speed  of  display  of  wave- 
forms on  our  storage  display  terminals.  We  are  currently  retrieving  SRO  data  from  the  Data- 
computer,  and  software  developments  in  this  area  are  described.  We  are  reducing  the  size  of 
our  library  collection  of  digital  seismic  data  so  that  proper  tape  maintenance  can  be  instituted. 

M.  A.  Chinnery 
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I.  ESTIMATION  OF  EVENT  LOCATION  AND  FOCAL  DEPTH 


A.  AUTOMATIC  EVENT  LOCATION  USING  3-COMPONENT 
SURF  ACE- WAVE  DATA 


The  association  of  surface-wave  trains  with  a specific  event,  particularly  when  successive 
surface-wave  trains  overlap,  is  a complicated  process.  Without  a reasonably  accurate  knowl- 
edge of  the  variation  of  group  velocity  across  the  surface  of  the  earth,  such  associations  are 
often  ambiguous  for  mixed  events.  We  describe  here  a technique  for  locating  seismic  events 
using  3-component  surface-wave  data,  which  may  also  be  of  use  in  locating  events  when  short- 
period  (SP)  data  are  not  available  because  of  very  high  levels  of  SP  noise  such  as  often  are  ob- 
served at  island  stations.  The  method  is  designed  to  be  used  in  conjunction  with  a power-based 
surface-wave  detector,  and  can  provide  preliminary  locations  with  a moderate  degree  of  accu- 
racy from  surface-wave  segments  as  short  as  60  sec  in  duration. 

Suppose  we  have  3-component  surface-wave  data  in  the  form  of  vertical,  north-south,  and 
east-west  seismometer  orientation,  denoted  by  Z(t),  N(t),  and  E(t),  respectively.  We  apply  a 
bandpass  filter  at  20  sec  period  to  each  channel,  and  use  a power  detector  to  determine  arrivals 
of  Rayleigh-wave  Airy  phases  on  the  vertical  component.  If  the  arrival  is  coming  from  an  azi- 
muth e with  respect  to  north,  then  radial  R(t)  and  transverse  Tr(t)  components,  given  by 


R(t)  = -cose  N(t)  - sin©  E(t) 


Tr(t)  = +sine  N(t)  - cose  E(t), 


will  then  contain  radial-component  Rayleigh  waves  and  transverse -component  Love  waves 
respectively. 

In  addition,  R(t)  and  Z(t)  are  polarized  such  that 
Z(t)  = ieR(t). . . 


where  e,  the  ellipticity,  is  a constant  for  a particular  frequency,  depending  only  upon  structure 
beneath  the  recording  station,  and  the  factor  i accounts  for  the  90°  phase  difference  between 
radial  and  transverse  components.  If  we  phase  shift  Z(t)  by  90°,  it  will  be  in  phase  with  R(t). 
Suppose  we  form  the  cross-correlogram 


of  R(t)  with  the  (90°  phase  shifted)  Hilbert  transform  Zn(t)  of  Z(t)  over  a time  window  of  length 
T,  then 


sine 


sine  I 


will  have  a maximum  value  at  the  epicenter  azimuth  e,  or  at 


given  by  © - tan 
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We  may  thus  determine  azimuth  6 via  Eq.  (1-4)  by  extracting  a window  of  length  T of  Z(t) 
centered  around  the  time  given  by  the  detector  output,  calculating  Z^(t).  and  cross-correlating 
it  with  the  same  windows  in  time  of  N(t)  and  E(t).  We  may  generate  radial  and  transverse  com- 
ponents R(t)  and  Tr(t)  corresponding  to  Z(t)  through  Eqs.  (I— 1 ) and  (1-2).  A check  upon  the  ac- 
curacy of  6 is  provided  by  the  ratio  of  R(t)/Z^(t)  at  peak  power  of  Z(t):  this,  through  Eq.  (1-3), 
should  be  the  ellipticity  *,  which  must  be  greater  than  — 1.0.  Typical  values  of  t at  20  sec 
period  are  -0.6  to  -0.9.  Now  that  we  have  Tr(t),  we  can  attempt  to  detect  the  Love-wave  ar- 
rival corresponding  to  the  Rayleigh  wave  detected  by  searching  for  maximum  power  of  Tr(t) 
over  a window  preceding  the  Rayleigh  wave.  The  Rayleigh-  and  Love-wave  arrival  times  T^ 
and  T^  are  given  by 

tr=4  and  Tl=^l 

where  A is  the  epicentral  distance,  and  U ^ and  are  the  Rayleigh-  and  Love-wave  group  ve- 
locities at  20  sec  period.  Now, 


TL  - TR  = A (u^  - U^)  = A(SL  - SR>  = 6AS 


where  S is  the  slowness  l/U.  While  both  and  vary  considerably  over  different  tectonic 
structures,  their  slowness  difference  6S  is  remarkably  constant;  in  a continuation  of  a study  of 
Rayleigh-  and  Love-wave  group  velocities  across  Eurasia,1  6S  values  of  0.037  with  a standard 
deviation  of  0.004  have  been  obtained  from  a suite  of  several  hundred  events. 

We  have  tested  this  technique  on  a suite  of  3-component  seismograms  recorded  at  Mashhad 
SRO  (MAIO)  from  78  events  within  Eurasia  at  distances  of  10“  to  75*.  The  propagation  paths  in- 
volved are  of  extreme  structural  complexity.  Tests  of  various  window  lengths  T for  the  deter- 
mination of  azimuth  0 showed  that  windows  as  short  as  60  sec  in  length  provided  the  best  results. 
We  used  6S  = 0.037  sec/km  for  the  distance  determination.  Figures  I-l(a)  and  (b)  show  the  er- 
rors in  azimuth  0 and  distance  A (deviations  from  actual  values  based  on  PDE  event  locations). 
Also  shown  in  Fig.  1-1  is  the  mislocation  (c),  defined  as  the  distance  between  the  true  epicenter 
and  that  of  a point  at  azimuth  0 and  distance  A from  the  station.  The  azimuth  is  determined 
to  within  15°  for  56  percent  of  the  79  events  studied;  but  an  alarmingly  high  27  percent  of  the 
events  are  mislocated  in  azimuth  by  more  than  30 “.  However,  if  we  apply  the  constraint  that 
the  ellipticity  e at  peak  power  lies  between  —0.5  and  —1.0  [satisfied  for  42  (or  53  percent)  of  the 
events  used],  this  reduces  the  number  of  such  gross  mislocations  dramatically.  The  distance 
determination  scheme  works  reasonably  well  (34  percent  within  5“),  and  60  percent  within  15“ 
of  true  distance;  here  again,  the  number  of  events  for  which  6A  exceeds  30“  is  reduced  by  ap- 
plication of  the  ellipticity  constraint  on  azimuth.  The  mislocation  statistics  are  not  very  encour- 
aging - only  50  percent  of  the  events  are  located  to  better  than  15“  in  distance.  These  large 
mislocations  are  primarily  due  to  errors  in  the  azimuth  determination,  and  an  improvement  of 
the  latter  is  clearly  required.  The  technique  requires,  however,  a minimal  amount  of  computer 
time  and  is  designed  to  he  used  in  conjunction  with  a surface-wave  detection  algorithm.  It  may 


prove  to  be  of  use  in  a surface-wave  association  scheme,  and  we  propose  further  investigation 
of  its  application.  R<G  North 

B.  THE  MASTER-EVENT  METHOD  FOR  LOCATING  EXPLOSION 

EPICENTERS  USING  CRUSTAL  PHASES 

This  report  concludes  the  experiment  of  evaluating  the  utility  of  crustal  phases  in  locating 

explosion  epicenters.  Previously,  we  determined  a two-layer  crustal  model  for  the  Basin  and 
2 3 

Range  province  ' within  which  the  majority  of  our  near-regional  stations  are  located.  Next, 
we  relocated  events  Greeley,  Scroll,  and  Rulison  using  the  Jeffreys-Bullen  (J-B)  travel-time 
tables  for  teleseismic  arrivals  and  our  own  Basin  and  Range  tables  for  crustal-phase  arrivals. 
These  relocation  results  showed  that  a large  amount  of  data  in  the  regional  distance  range  was 
able  to  produce  epicenters  as  accurate  as  those  derived  from  substantial  teleseismic  data.  Fur- 
thermore, adding  observations  from  a few  good  regional  stations  to  a teleseismic  data  set  usually 
improves  the  resulting  epicenter  in  both  accuracy  and  quality.  Finally,  to  estimate  the  effect 
of  bias  in  the  travel-time  tables,  we  present  here  the  results  of  relocating  the  events  Rex,  Box- 
car, and  Commodore  using  the  master-event  method.  This  method  is,  in  effect,  a means  of 
automatically  applying  the  usual  station  corrections. 

The  NTS  explosion  Greeley  was  selected  as  the  master  event  because  it  had  the  best  distri- 
bution of  reporting  stations.  The  corresponding  travel-time  residuals  were  then  subtracted 
from  the  arrival  times  for  the  other  events.  The  difference  between  the  Greeley  residual  and 
the  average  residual  of  all  four  explosions  averaged  over  the  set  of  common  stations  had  a mean 
of  0.68  sec,  with  a standard  deviation  of  4.9  sec.  This  shows  that  Greeley  residuals  are  repre- 
sentative of  all  four  explosions. 

The  results  of  relocating  the  epicenters  of  Rex,  Boxcar,  and  Commodore  after  correcting 
for  the  Greeley  residuals  are  listed  in  Tables  1-1,  1-2,  and  1-3.  They  show  that  the  master-event 
method  gives  more-accurate  epicenters  than  by  using  uncorrected  arrivals  in  all  three  cases. 

The  reduction  in  mislocation  varies  from  approximately  a factor-of-2  to  a factor-of-6  and  de- 
pends on  the  number  of  reporting  stations.  However,  there  does  not  appear  to  be  any  significant 
improvement  in  the  quality  of  the  epicenter  as  measured  by  the  area  of  the  corresponding  confi- 
dence ellipse.  In  the  case  of  a large  well-recorded  event  such  as  Boxcar,  these  results  show 

that  the  epicenter  can  be  located  approximately  1 km  away  from  its  known  position  with  a confi- 

2 

dence  ellipse  area  of  less  than  8 km  . 

We  wish  to  emphasize  a general  conclusion  of  these  location  experiments.  Accurate  epi- 
centers (<5  km  mislocation)  can  be  derived  from  crustal-phase  data  if  either  station  corrections 
are  known  or  if  the  crustal  structure  is  both  simple  and  well  determined.  Of  the  two  options, 
station  corrections  measured  by  observing  a master  event  are  preferable  to  just  knowing  the 
crustal  structure.  However,  in  an  area  as  homogeneous  and  well  known  as  the  Basin  and  Range, 
epicenters  accurate  to  within  10  km  can  be  found  without  resorting  to  a master  event.  The  ad- 
vantage to  the  latter  approach  is  that  the  epicenters  will  stay  accurate  over  a range  of  source 
and  station  locations  as  long  as  they  both  remain  within  the  same  geologic  province.  Strictly 
speaking,  station  corrections  measured  from  a master  event  will  only  apply  to  events  in  the 
immediate  vicinity  of  that  epicenter.  R>  E>  Needham 

D.  W.  McCowan 
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C.  SEISMIC  RAY  TRACING  FOR  RELOCATING  OCEANIC  RIDGE  EARTHQUAKES 


The  wide  scatter  and  systematic  offset  of  earthquake  epicenters  away  from  the  narrow  tec- 
tonic axis  defined  by  the  oceanic  ridge  crests  suggests  that  large  location  errors  may  exist. 

The  errors  could  result  from  the  fact  that  oceanic  ridge  earthquake  events  are  small  (mb  < 5.0) 
and  thus  are  only  recorded  by  nearby  stations  whose  azimuthal  distribution  is  generally  biased. 
Accordingly,  for  earth- structure  models  which  have  an  average  velocity  slower  than  the  true 
velocity,  the  location  will  be  "pulled"  along  the  general  azimuth  toward  the  region  of  the  greatest 
number  of  stations.  Conversely,  a "fast”  model  would  "push"  the  location  away  from  the  domi- 
nant stations'  region.  Also,  azimuthal  ray  bending  due  to  the  shallow  two-dimensional  structure 

beneath  the  ridges  could  introduce  errors. 

2 3 

In  previous  SATS,  ’ Phillips  and  others  have  employed  realistic  earth  velocity  models  of 
oceanic  ridge  structure  to  calculate  travel-time/epicentral-distance  tables.  Ray-tracing 

4 

methods  were  used  to  account  for  both  the  higher  crustal  velocity  and  two-dimensional  earth 
structure  beneath  oceanic  ridges.5  Significantly,  the  new  ocean  ridge  model  produced  much 
faster  travel  times  as  compared  with  the  spherical  earth,  Jeffrey-Bullen,^  and  Herrin7  models 
which  are  generally  used  to  locate  earthquakes  [i.e..  International  Seismological  Center  (ISC)). 

As  might  be  expected,  the  most  marked  travel-time  difference  occurred  for  near-surface 
events.  Travel  times  were  typically  4 sec  faster  than  those  computed  with  the  J-B  model  over 
the  epicentral  distance  range  20°  to  50°.  Between  50°  and  70°  the  difference  was  5 sec,  and  for 
distances  less  than  20°  the  ocean  ridge  model  computed  travel  times  were  6 sec  faster.  Notably, 
the  dependence  of  travel  time  on  the  azimuth  to  the  station  relative  to  the  ridge's  trend  is  very 

g 

small  (<0.3  sec)  for  epicentral  distances  greater  than  10°.  Solomon  and  Julian  drew  a similar 
conclusion  in  their  study  of  azimuthal  ray-bending  effects  on  focal-plane  solutions  along  ridges. 

The  effect  of  using  a realistic,  oceanic  ridge  velocity  model  for  relocating  mid-Atlantic 
ridge  earthquakes  was  demonstrated  in  a previous  SATS3  for  a few  events  (10)  in  a small  region 

north  of  the  Azores.  Here,  the  ridge's  tectonic  axis  is  especially  well  defined  by  high-precision 
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bathymetric  and  micro-earthquake  surveys.  This  work  clearly  showed  that  the  use  of  travel- 
time  tables  derived  from  the  faster  ocean  ridge  earth  model  resulted  in  moving  the  ISC-reported 
locations  westward,  away  from  the  dominant  North  American  station  region  and  much  closer  to 
the  ridge' s tectonic  axes.  In  fact,  by  depth  constraining  the  events  to  the  base  of  the  oceanic 
lithosphere  (~110  km),  most  events  could  be  moved  exactly  onto  the  ridge  crest  axis!  However, 
since  the  maximum  depth  of  events  along  the  spreading  axis  is  probably  less  than  10  km  due  to 
the  thin  lithosphere  here,  such  a depth  constraint  is  not  reasonable.  Alternatively,  the  same 
relocation  effect  could  also  be  accomplished  by  increasing  the  average  velocity  of  the  ocean 
ridge  model  even  more  over  the  J-B  model  to  provide  faster  travel  times.  An  additional  travel- 
time decrease  of  about  1 sec  (25  percent)  for  the  epicentral  distance  range  of  20°  to  70°  for  a 
near-surface  event  would  approximate  the  110 -km  depth  constraint. 

We  have  recently  completed  a broader  test  of  the  efficacy  of  the  ocean  ridge  velocity  model 
(ORM)  for  relocating  earthquakes.  For  this  work,  we  chose  two  large  regions  on  the  North  and 
Central  Atlantic  ridges  (Figs.  1-2  and  1-3)  where  there  is  pronounced  systematic  offset  and  wide 
scatter  of  the  ISC  earthquake  epicenters  relative  to  the  well-charted  ridge  axes;  86  events  were 
examined  in  all  {42  and  44,  respectively).  The  azimuth  to  the  dominant  recording-station  region 
in  North  America  for  both  areas  is  to  the  northwest.  The  relocation  program  of  the  DADS  for 
the  PDP-7  computers  was  used  to  individually  process  each  event. 
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TABLE  1-4 


* 


The  results  of  these  relocations  are  shown  by  the  arrows  in  Figs.  1-2  and  1-3.  Note  the  re- 
duced scatter  and  general  shift  of  the  earthquake  epicenter  pattern  toward  the  ridge  axis.  These 
effects  are  summarized  quantitatively  in  Table  1-4  where  the  mean  values  of  the  offset  distance 
for  the  J-B  and  ORM  models  are  compared.  The  decrease  for  both  regions  combined  is  about 
4 km  (i.e.,  14.1  vs  10.2  km)  for  the  ORM.  For  the  Central  Atlantic  epicenter,  it  is  even  more 
dramatic  (i.e.,  15.52  vs  9.8  km).  Also,  the  change  in  the  variance  of  the  mean  offset  values  is 
significant  at  the  9 5-percent  probability  for  both  regions. 

The  eastward  and  southward  shift  of  each  region's  epicenter  pattern,  respectively,  is 
demonstrated  by  the  decrease  in  the  accumulative  offset-distance  values  in  Table  1-4.  That  is, 
if  seismicity  is  normally  distributed  about  the  plate  boundary  defined  by  the  ridge  axis,  as  might 
be  expected,  the  sum  of  the  offset  distances  should  be  zero  for  a large  number  of  events.  A 
large  accumulative  distance  value  would  suggest  a systematic  offset.  Note  that  for  both  the  J-B 
and  ORM  models,  the  locations  are  to  the  west  and  north  of  the  plate  boundary.  However,  in 
both  regions  the  accumulative  offset-distance  values  decrease  markedly  toward  zero  for  the 
ORM  relocated  events.  In  fact,  there  has  been  an  approximately  70-percent  decrease  in  the 
accumulated  offset  distance  (i.e.,  248  vs  72  km  and  95  vs  28  km,  respectively).  Presumably, 
an  additional  30 -percent  faster  travel-time  decrease  for  the  ORM  over  the  J-B  travel  times 
would  have  centered  the  epicenter  pattern  over  the  ridge  axis.  Such  a model  would  imply  an 
additional  l-sec  decrease  in  travel  times  as  compared  with  the  J-B  model  derived  tables.  This 
is  in  reasonable  agreement  with  the  110  km  depth  constrained  relocation  effect  noted  earlier. 

As  a result  of  these  encouraging  test  results,  we  are  now  developing  a relocation  program 
for  the  PDP-11  computers  to  batch  process  all  mid-Atlantic  ridge  earthquake  epicenters  with 
travel-time  tables  derived  from  the  new  ocean  ridge  structure  model.  This  work  should  provide 
a much-improved  delineation  of  the  ridge  axis,  especially  in  regions  such  as  the  South  Atlantic, 
where  little  topographic  information  is  presently  available. 

J.  D.  Phillips 
R.  E.  Needham 
R.  M.  Sheppard 

D.  REI.ATIVE  DEPTH  AND  SEISMIC  MONITORING 

It  seems  a reasonable  proposition  that  all  major  shallow  earthquakes  occur  below  the  depth 
of  the  deepest  underground  explosions.  If  this  were  true,  then  earthquake  monitoring  in  a given 
region  would  be  simplified  b;  accurate  location  of  each  event  in  a region  with  respect  to  the 
shallowest  large  earthquake  in  the  same  region.  Only  those  events  significantly  above  the  mas- 
ter earthquake  would  require  special  attention.  Such  a culling  scheme  is  only  practical  in  a 
tectonically  active  region  such  as  Central  Asia  or  along  the  Kuril  and  Kamchatka  arcs.  How- 
ever, it  is  in  these  regions  where  culling  is  essential  to  make  the  monitoring  effort  manageable, 
given  the  immense  quantity  of  incoming  data. 

To  gain  some  insight  into  the  practical  difficulties  of  using  such  a scheme,  two  data  sets 
were  compiled  from  ISC  arrival  times  stored  in  a computer  disk  file.  One  set  was  from  after- 
shocks of  an  earthquake  that  occurred  beneath  the  Aleutian  trench  (Fig.  1-4)  about  2 months  after 
the  Rat  Island  earthquake  of  4 February  1965.  The  other  set  was  from  earthquakes  in  the  south- 
ern part  of  the  Kuril  arc  (Fig.  1-5),  one  of  the  world's  most  active  regions  of  shallow  earthquakes 
in  the  last  15  years. 
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The  relative  locations  were  computed  by  a technique  described  in  Fitch  and  Jackson,  and 
Jackson  and  Fitch.2  The  uncertainties  in  relative  depth  are  primarily  controlled  by  the  station 
coverage  for  each  secondary  event.  Nonrandom  errors  in  reading  arrival  times  can  signifi- 
cantly bias  the  relative  location  and,  in  particular,  the  relative  depth  of  the  less-well-recorded 
events.  To  some  extent,  this  is  a consequence  of  using  a least-squares  of  L.2  norm  as  a solu- 
tion criterion. 

The  Rat  Island  locations,  shown  in  vertical  sections  normal  to  the  arc  in  Figs.  1-6  and  1-7, 
reveal  activity  in  a zone  inclined  toward  the  north.  The  station  distribution  is  grossly  non- 
uniform,  with  the  less-well-recorded  aftershocks  having  no  station  coverage  in  the  southern 
half  of  the  focal  sphere.  Without  arrival  times  from  ADAK,  the  closest  station  at  a distance  of 
about  3 50  km  toward  the  ENE,  these  locations  are  significantly  more  scattered  about  the  master. 
The  importance  of  ADAK  is  illustrated  in  Figs.  1-6  and  1-7.  A much  shallower  inclination  to  the 
zone  and  diminished  uncertainties  result  from  a change  from  56°  to  9 5°  in  the  takeoff  angle  for 
the  ray  path  from  the  master  to  ADAK.  This  is  equivalent  to  assuming  a direct  rather  than  a 
refracted  first  arrival  at  ADAK.  In  either  case,  about  half  the  activity  is  above  the  master  if 
the  measures  of  uncertainty  are  taken  as  one  standard  deviation  which  is  more  than  20  km  for 
the  less-well-recorded  aftershocks.  To  be  useful  in  the  context  of  seismic  monitoring,  smaller 
uncertainties  are  required  — say,  less  than  10  km  for  one  standard  deviation. 

Figure  1-8  shows  that  the  majority  of  the  Kuril  earthquakes  are  better  located  with  respect 
to  their  master  than  are  the  Rat  Island  aftershocks.  Station  coverage  is  more  complete  because 
of  the  large  number  of  Japanese  stations  and  a significant  increase  in  the  total  number  of  report- 
ing stations  in  the  5-year  interval  between  the  times  of  the  Rat  Island  and  Kuril  master  earth- 
quakes. In  fact,  one  standard  deviation  in  relative  depth  is  less  than  10  km  for  Kuril  events 
located  with  more  than  60  differential  arrival  times.  This  is  illustrated  in  Fig.  1-9.  In  Fig.  1-8, 
the  events  with  the  larger  uncertainties  form  a halo  on  the  seaward  side  of  the  dense  cluster  of 
events  that  mark  the  contact  zone  between  the  Asian  and  Pacific  plates.  Apparently,  the  events 
in  the  halo  have  biased  as  well  as  highly  uncertain  locations  relative  to  the  master.  These  events 
account  for  15  percent  of  the  total  number. 

If  station  coverage  improves  with  time,  as  it  surely  has  in  the  8 years  since  the  Kuril  mas- 
ter earthquake,  the  utility  of  a relative  location  scheme  such  as  the  one  proposed  here  must  also 
improve.  The  same  cannot  be  said  for  earthquake  locations  computed  with  uncalibrated  travel 

tUnes-  T.J.  Fitch 
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Fig.  1-2.  Chart  showing  relocation  (solid  circles)  of  ISC  Bulletin  epicenters 
(open  circles)  for  North  America  and  Europe  using  ocean  ridge  structure  ve- 
locity model  (ORM).  Note  closer  grouping  of  epicenters  about  ridge  axis  for 
ORM  relocated  events.  Events  are  for  period  1964-1973.  Solid  line  shows 
ridge  axis. 
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Fig.  1-6.  Rat  Island  relative  location:  56°  takeoff  angle  for  ray 
path  to  ADAK. 
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Fig.  1-8.  Kuril  relative  locations:  uncertainties  on  one  standard  deviation. 


Fig.  1-9.  Size  of  Kuril  data  sets  vs 
uncertainty  in  relative  depth. 
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II.  ESTIMATION  OF  SOURCE  MOMENT  TENSOR 
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A.  INVERSION  OF  BODY-PHASE  AMPLITUDES  FOR  THE  MOMENT  TENSOR 

A mid-Atlantic  earthquake  that  occurred  on  23  October  1964  was  chosen  for  this  inversion 
experiment.  The  source  region  is  beneath  an  abyssal  plain  and  therefore  free  of  major  struc- 
tural inhomogeneities.  The  body  waves  P,  pP.  and  sP  were  well-recorded  on  the  LP  instru- 
ments of  the  WWSSN  at  a time  when  their  seismometers  had  a 30-sec  rather  than  the  present 
15-sec  free  period.  The  azimuthal  coverage  afforded  by  stations  in  Europe,  North  America, 
and  Africa  is  as  complete  as  can  be  expected  from  a teleseismic  data  set. 

The  analysis  of  the  analog  records  proceeded  as  follows:  I.P  vertical  seismograms  from 
27  stations  in  the  teleseismic  distance  range  from  25°  to  90°  were  digitized.  DC  and  linear 
biases  were  removed,  and  linear  interpolation  was  used  to  resample  the  data  at  a constant  time 
increment  of  0.25  sec.  Examples  of  the  resampled  data  are  shown  in  Fig.II-1.  To  measure  the 
relative  excitations  of  p,  pP,  and  sP  the  digitized  seismograms  were  spiked  by  convolution  with 
an  inverse  filter  that  removes  the  effect  of  the  instrument  and  the  anelasticity  of  the  earth  equiv- 
alent to  a t*  of  0.5.  The  effect  of  the  inverse  filter  is  illustrated  in  Figs.II-2  through  II-4.  The 
construction  of  filters  of  this  kind  is  discussed  in  Robinson.1  The  raw  filtered  records  contained 
high-frequency  noise  that  is  thought  to  be  an  unavoidable  consequence  of  digitizing.  To  minimize 
this  noise,  various  low-pass  filters  were  applied  to  the  spiked  records.  The  most  satisfactory 
one  was  a phaseless  filter  with  a 12-dB-per-octave  falloff  at  1 Hz.  Examples  of  these  records 
are  shown  in  Fig.II-5  with  the  phases  identified.  ' 

The  data  for  the  moment  tensor  inversion  were  the  amplitude  of  the  appropriate  spikes. 

The  mean  pP-P  time  of  10.5  sec  and  sP-P  time  of  14.7  sec  are  consistent  with  a source  depth 
of  about  45  km  below  the  seafloor.  In  Fig.  II-6  it  is  apparent  that  there  is  excellent  coverage  of 
one  compressional  quadrant  of  the  focal  sphere.  The  data  set,  80  amplitudes,  is  sufficient  to 
resolve  significant  differences  from  a standard  double-couple  source  if  such  differences  exist. 

For  moment  tensor  inversion,  the  amplitudes  must  be  normalized  for  distance  and  free- 
surface  effects.  The  distance  effects  include  geometric  spreading  and  anelasticity.  The  latter 
has  been  removed  in  an  ad  hoc  manner  by  including  a t*  in  the  construction  of  the  inverse  filter. 
The  former  has  been  removed  by  first  normalizing  all  amplitudes  to  an  arc  distance  of  60°  with 
an  amplitude  distance  function  derived  from  LP  data  presented  in  Sengupta  and  Julian.^  A stan- 
dard spreading  factor  for  a distance  of  60°  was  then  used  to  reduce  amplitude  to  the  focal 
sphere.3  The  near-source  conversion  of  s to  P and  p to  P and  the  near-station  effect  of 
the  free  surface  were  removed  oy  applying  geometric  ray  theory  to  the  case  of  plane  waves 
incident  on  a free  surface  from  a half  space.  The  details  of  this  part  of  this  analysis  follow 
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Kanamori  and  Stewart. 

The  reduced  amplitudes  were  inverted  for  the  moment  tensor  source  representation  using 
the  formula  derived  by  McCowan.  The  solution  revealed  a source  that  contained  90-percent 
double  couple,  9-percent  compensated  linear  vector  dipole,  and  1 -percent  monopole  — essen- 
tially a pure  double  couple  as  is  assumed  in  first-motion  solutions.  The  reduced  amplitudes 
fit  the  model  with  an  rms  of  16  percent.  However,  in  this  case  the  nodal  planes  are  oriented 
about  90°  counterclockwise  from  planes  inferred  from  a first-motion  solution  by  Sykes  and 
Sbar.^  The  contrasting  double  couples  are  summarized  in  Table  II- 1 and  illustrated  in  Fig.  II-6. 
The  nodal  amplitudes  recorded  at  the  three  closest  stations  — BEC,  SJG,  and  TRN  — are,  at  a 
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TABLE  ll-l 

DOUBLE  COUPLE  FROM  MOMENT  TENSOR  INVERSION 


Strike 

Slip 

Dip 

Nodal  Plane  1 

309 

292 

50 

Nodal  Plane  2 

160 

246 

44 

Trend 

Plunge 

Pole  of  Plane  2 

70 

45 

Pole  of  Plane  1 

219 

40 

P-axis 

53 

4 

T-axis 

152 

73 

Sykes  and  Sbar  (Ref. 6) 

Strike 

Slip 

Dip 

Nodal  Plane  1 

44 

315 

57 

Nodal  Plane  2 

286 

221 

53 

Trend 

Plunge 

Pole  of  Plane  1 

314 

33 

Pole  of  Plane  2 

195 

37 

P-axis 

164 

3 

T-axis 

258 

54 

first  glance,  more  consistent  with  the  Sykes  and  Sbar  solution;  but  the  uncertainties  in  the 
takeoff  angles  to  these  stations,  all  in  the  distance  range  9°  to  15°,  are  sufficiently  great  that 
consistency  with  the  inversion  double  couple  can  not  be  precluded.  The  P-axes  (representing 
compressive  stress)  for  both  mechanisms  are  nearly  horizontal;  however,  the  orientation  of 
the  P-axis  from  the  inversion  double  couple  is  more  nearly  normal  to  the  mid-Atlantic  ridge. 

Such  an  orientation  is  preferred  by  focal  mechanisms  of  the  thrust  type  in  eastern  North  America 
and  beneath  the  North  Atlantic  basin.6  Inversion  for  the  full  moment  tensor,  as  was  done  in  this 
experiment,  should  yield  source  components  that  are  less  influenced  by  systematic  errors  in  the 
reduced  amplitudes  than  inversions  for  a constrained  moment  tensor  as  proposed  by  Strelltz.7 

T.J.  Fitch 
M.  W.  Shields 

t D.  W.  McCowan 

B.  STRUCTURAL  INHOMOGENEITY  AND  THE  RESOLUTION 

OF  SOURCE  FUNCTION 

. 

In  a recent  SATS,5  Frasier  and  Fitch  showed  that  long-period  (LP)  P-amplitudes  from 
Cannikin  varied  as  much  from  station- to- station  as  did  the  short-period  (SP)  amplitudes.  The 
simplest,  albeit  indeterminate,  explanation  for  this  result  Is  transmission  of  the  seismic  waves 
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through  and  around  a structurally  inhomogeneous  source  region  for  which  there  is  ample  geo- 
physical evidence  Here,  we  investigate  the  effect  of  such  amplitude  distortion  on  inversions 
for  the  moment  tensor  source. 

The  Cannikin  LP  amplitudes  in  Fig.II-7,  which  are  normalized  to  a distance  of  60°,  are 
consistent  with  a moment  tensor  source  with  36-percent  monopole,  57-percent  double  couple, 
and  7 -percent  compensated  linear  vector  dipole  (CLVD).  The  double-couple  component  is 
superposed  on  the  amplitude  distribution  in  Fig.  II-7.  The  degree  of  amplitude  distortion  from 
near-source  structure  is  probably  near  an  upper  limit  for  this  case.  The  frequency  content  of 
these  signals  is  concentrated  near  the  high-frequency  end  of  the  teleseismic  source  spectrum 
in  the  bandpass  1 to  0.2  Hz,  and  the  structural  inhomogeneity  associated  with  the  Aleutian  arc 
and  subduction  may  involve  velocity  contrasts  of  about  5 to  1 0 percent  over  horizontal  distances 
of  tens  to  hundreds  of  kilometers. 

To  add  credibility  to  the  interpretation  of  the  Cannikin  amplitudes,  a similar  analysis  was 
carried  out  with  P-wave  amplitudes  from  an  intermediate  earthquake  located  about  5°  east  of 
the  test  site  on  Amchitka.  Figure  II-8  shows  that,  as  in  the  case  of  the  Cannikin  amplitudes, 
those  for  the  earthquake  show  nearly  as  much  LP  or  SP  variation.  In  this  case,  the  LP  am- 
plitudes are  dominated  by  spectral  components  near  the  15-sec  peak  in  instrument  response. 

The  apparent  moment  tensor  representation  is  60-percent  double  couple,  29-percent  CLVD. 
and  11 -percent  monopole.  A large  compressional  amplitude  at  College  Station  (COL  in  Fig.  II— 9 ) 
accounts  for  about  one-third  of  the  CLVD  component.  Instrument  polarity  at  this  station  has 
been  found  to  be  correct.  The  non-double-couple  components  are  likely  to  be  artifacts  of  an 
incomplete  data  set  and  amplitude  distortion.  Figure  II-9  shows  the  distribution  of  first-motion 
and  the  double-couple  component. 

If  the  same  analysis  is  applied  to  LP  P-amplitudes  from  sources  in  structurally  less- 
complex  parts  of  the  world,  the  double  couple  (in  the  case  of  earthquakes)  and  the  monopole 
(in  the  case  of  explosions)  dominates  the  moment  tensor  representations.  For  example,  the 
Borrego  Mountain  earthquake  in  Southern  California  yielded  81-percent  double  couple,  13-percent 
CLVD,  and  6-percent  monopole.  The  plot  of  LP  vs  SP  amplitudes  in  Fig.  II- 1 0 shows  a distinct 
elongation  along  the  SP  axis  which  is  expected  in  cases  where  large-scale  structural  complexity 
is  absent  in  the  same  region.  The  amplitude  distribution  and  double-couple  component  are  shown 
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by  McCowan  elsewhere  in  this  SATS.  A Novaya  Zemlya  explosion  has  a much-more-uniform 
amplitude  distribution  than  Cannikin  (Fig.  II— 1 1 ) . The  moment  tensor  components  are  62-percent 
monopole,  34-percent  double  couple,  and  3-percent  CLVD. 

In  this  study,  loose  upper  bounds  have  been  placed  on  the  effect  of  amplitude  distortion  on 
moment  tensor  sources.  Structurally  related  distortion  can  not  be  removed  by  any  known 
scheme.  Consequently,  such  distortion  remains  a major  obstacle  to  the  resolution  of  source 
functions,  particularly  where  sources  are  located  in  structurally  complex  regions. 

T.J.  Fitch 
M.  W.  Shields 

C.  ESTIMATING  FAULT-PLANE  SOLUTIONS  BY  ROBUST  METHODS 

Finding  fault-plane  solutions  for  earthquakes  (or  explosions)  remains  one  of  the  last  ex- 
amples of  manual  drudgery  in  seismological  research.  Unfortunately,  the  procedure  is  a 
familiar  one.  First  the  data,  usually  LP  polarities,  are  read  either  from  film  or  digital 


recordings.  Next,  each  data  point  is  plotted  on  a stereographic  projection  of  the  lower  focal 
hemisphere.  Then,  the  fault-plane  solution  is  found  by  manipulating  a stereonet  on  top  of  the 
data  projection.  Usually,  there  is  a spectrum  of  allowable  fault-plane  solutions  which  will 
satisfy  the  data  so,  finally,  bounds  on  the  solution  are  estimated  by  mapping  out  extreme  values 
with  the  stereonet. 

The  objective  of  this  research  is  twofold:  first,  to  reduce  the  considerable  labor  involved; 
and  second,  to  remove  (as  much  as  possible)  the  subjective  element.  Both  these  difficulties 
are  common  sources  of  error.  It  is  hoped  that  removing  them  will  improve  routine  fault-plane 
solutions  and  therefore  lead  to  better  discrimination  between  earthquakes  and  explosions. 

The  operational  problems  involved  in  finding  fault-plane  solutions  are  threefold.  First, 
the  data  consist  only  of  inequalities  which  will,  of  course,  admit  a wide  variety  of  models. 
Second,  the  data  must  be  corrected  by  projecting  them  back  through  the  transmitting  earth  to 
the  focal  sphere,  which,  of  necessity,  implies  the  assumption  of  an  earth  model.  Third,  the 
data  are  satisfied  by  an  ad  hoc  trial-and-error  procedure  of  fitting  the  model,  which  involves 
an  element  of  subjective  judgment. 

The  first  problem  in  the  procedure  can  be  solved  easily  by  using  the  LP  P-wave  amplitudes 
themselves  instead  of  just  their  signs.  There  is,  it  turns  out,  very  little  extra  labor  involved 
in  measuring  and  recording  the  relative  amplitudes  on  paper  as  there  is  in  noting  their  signs. 
Unfortunately,  the  second  problem  then  becomes  significant.  Finding  where  a ray  took  off  from 
the  focal  sphere  to  a given  station  is  far  easier  than  correcting  its  amplitude  observed  on  the 
surface  of  the  earth  back  to  the  source.  On  the  other  hand,  if  properly  corrected  amplitudes 
are  available,  then  the  third  problem  is  virtually  solved  because  deterministic  estimation 
schemes  can  be  readily  employed. 

So,  the  major  stumbling  block  would  appear  to  be  in  correcting  the  amplitude  observations 
back  to  the  source.  In  fact,  there  is  another  small  problem  which  places  a constraint  on  any 
practical  method.  Frequently,  due  to  variations  in  signal-to-noise  ratio  (S/N),  it  is  not  possible 
to  measure  an  amplitude  accurately  but  it  is  still  possible  to  observe  its  polarity.  Any  workable 
method  must  therefore  allow  data  that  can  be  either  amplitudes  or  just  signs. 

We  attempted  to  solve  the  second  problem  by  adopting  an  estimation  procedure  based  on 
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minimizing  the  Lt  norm.  The  method,  as  formulated  by  Claerbout  and  Muir,  is  a way  of  find- 
ing the  optimum  linear  model  which  minimizes  the  sum  of  the  absolute  differences  between  it- 
self and  the  data.  This  is  opposed  to  the  usual  least-squares,  or  L2,  norm  which  minimizes 
the  sum  of  the  squares  of  the  corresponding  differences.  The  LI  norm  has  two  important  ad- 
vantages over  the  L2  norm:  first,  it  is  less  sensitive  to  erroneous  data  values,  a property 

called  robustness;  and  second,  it  conveniently  allows  data  that  consist  only  of  inequalities.  We 
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applied  the  method  by  fitting  the  linear  point  source  moment  tensor  model  directly  to  the  film- 
chip  amplitudes  and  polarities.  In  this  way,  correcting  those  amplitudes  back  to  the  source 
was  completely  neglected.  Because  of  its  robust  properties,  the  LI  norm  estimates  linear 
models  which  are  less  sensitive  to  such  a gross  omission  as  would  be  those  estimated  by  min- 
imizing an  L2  norm. 

To  illustrate  the  method,  we  chose  two  events:  the  Borrego  Mountain  earthquake  of  9 April 
1968,  and  a deep  Fiji-Tonga  earthquake  of  9 October  1967.  An  analyst  read  the  LP  amplitudes 
for  both  events  from  WWSSN  film  chips.  Figures  11-12  and  11-13  show  equal-area  projections 
of  the  lower  focal  hemisphere  of  the  Borrego  Mountain  earthquake  with  the  data  and  fault  planes 
indicated  for  the  LI  and  L2  norms,  respectively.  In  the  latter  case  (as  in  the  other  L2  norm 
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examples  given  below),  the  solution  was  constrained  to  have  no  volume  change.  Both  solutions 
are  similar  in  that  they  show  a predominance  of  double-couple  source  mechanism  with  a strike 
slip  orientation.  Figures  11-14  and  11-15  show  the  same  cases  except  that  the  amplitude  of  one 
data  point  has  been  multiplied  by  ten.  As  can  be  seen,  the  LI  solution  is  still  much  the  same 
as  it  was,  but  the,  L2  solution  is  completely  different  — it  is  still  dominated  by  a strike  slip 
double  couple  but  has  been  rotated  by  nearly  45°.  This  example  illustrates  the  robust  property 
of  the  LI  estimate. 

Figures  11-16  through  11-18  show  results  for  the  Fiji-Tonga  earthquake.  Figure  II- 16  is  an 
LI  norm  estimate  using  all  the  LP  amplitude  and  polarity  data,  and  F'ig.  11-17  is  an  LZ  norm 
solution  using  just  the  LP  amplitudes.  Both  are  similar,  although  the  1.2  norm  solution  has 
significantly  more  double-couple  mechanism  in  it.  Figure  11-18  is  the  same  as  Fig.  II - 16,  but 
with  an  additional  nine  near-in  SP  polarities  added.  Clearly,  this  solution  is  much  better  in 
terms  of  resolution  than  the  other  two;  it  is  composed  of  principally  double-couple  source 
mechanism.  So,  the  inclusion  of  a few  additional  polarities  completely  changes  a fault-plane 
solution  both  in  terms  of  its  orientation  and  force  composition.  If  we  were  to  restrict  ourselves 
to  using  just  amplitude  data,  the  inclusion  of  these  critical  data  points  would  not  have  been 
possible. 

These  results  suggest  that  a balanced  approach  to  calculating  fault-plane  solutions  based 
on  minimizing  the  Li  norm  of  LP  amplitude  residuals  supplemented  by  reliable  first-motion 
(polarity)  data  may  be  practical  in  general  usage.  It  is  a convenient  way  of  blending  in  the 
first-motion  data,  and  is  relatively  insensitive  to  erroneous  data  values. 

D.  W.  McCowan 

D.  SYNTHETIC  SEISMOGRAMS  BY  SUPERPOSITION  OF  NORMAL  MODES: 

APPLICATION  TO  THE  INVERSE  PROBLEM  FOR  SOURCE  MECHANISM 

Current  expansion  of  the  world-wide  network  of  stations  with  digital  recording  should  make 
possible  routine  determination  of  the  moment  tensor  of  seismic  sources.  The  effects  of  lateral 
heterogeneities  may  impose  certain  restrictions  on  the  upper  limit  of  the  frequency  range.  Also, 
it  may  be  necessary  to  avoid,  at  first,  surface  waves  whose  dispersion  is  a strong  function  of 
the  structure  in  the  most  heterogeneous  part  of  the  earth. 

The  purpose  here  is  to  describe  an  algorithm  that  can  be  used  to  efficiently  generate  syn- 
thetic seismograms  by  superposition  of  normal  modes,  and  to  show  application  of  the  method 
to  the  direct  and  inverse  problem  of  importance  in  seismology.  Although  only  the  modes  that 
were  obtained  by  the  exact  solution  will  be  used  here,  the  method  can  be  easily  combined  with 
results  of  computations  at  higher  frequencies  by  approximate  methods  that  are  faster  and  suf- 
ficiently accurate. 

Q 

Following  Gilbert  and  Dziewonski,  the  displacement  function  in  epicentral  coordinates  can 
be  written  as 

6 oo  oo 

u(r,t)  = £ M^t)*  Z l nf/(t)  £ ns/m<r)  * ln?/m‘ro'l  1 . 

i=l  n=0  1=0  m 

where  M^(t)  are  the  six  independent  components  of  the  moment  tensor  as  a function  of  time  and 
nff(t)  ={l-cosnwft.  exp (-n <*>,1/2^1}  • H(t) 


describes  the  temporal  behavior  of  each  individual  mode.  Summation  over  m extends  from 
— min  (2,/ ) to  min  (2,1 ).  Our  notation  follows  that  in  Ref.  9,  and  will  not  be  fully  explained  here. 
We  note  that  the  effect  of  the  geographical  coordinates  of  the  receiver  is  the  same  for  all  the 
modes  of  the  same  angular  order  I.  Thus,  for  a given  mode  type  and  angular  order  number, 
contribution  of  all  the  overtones  can  be  evaluated  first  — the  corresponding  functions  will  depend 
on  the  source  radius  r , receiver  radius  (usually  equal  to  the  radius  of  the  earth  — a),  and  time. 
For  spheroidal  modes,  one  needs  to  evaluate  four  functions  for  vertical  as  well  as  horizontal 
displacements: 
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For  toroidal  modes  there  are  only  two  functions,  denoted  and  W4  because  of  their 
analogy  with  V3/  and  V4 : 

W3l  = 2nW,<a»  Wro>  “ nWi(ro>/rl  nfl(t) 
n 

W.,  = r*1  V W,  (a)  W,  (r  ) f,  (t)  . 

41  o iJ  n t'  ' n t o n t'  ' 

n 

Examples  of  functions  U^,  V^,  and  computed  for  a source  depth  of  650  km  and  1 = 4 are 

shown  in  Fig.  11-19.  The  appropriate  eigenfunctions  and  eigenfrequencies  were  computed  by 
10 

Buland  for  the  earth  model  1066B  (Ref.  9);  the  quality  factors  Q,  were  obtained  using  the 
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anelastic  structure  of  the  model  FSQMK  of  Dzlewonski.  * The  time  traces  represent  super- 
position of  the  modes  with  periods  from  1545  sec  (gS4l  to  45  sec  Because  of  the  low 

angular  order  number,  at  high  frequencies  the  compressional  waves  can  be  associated  with 
PKIKP  close  to  180*  and  shear  waves  with  nearly  vertical  ScS  reflections.  The  sharp  onset 
on  Uj  is  the  arrival  of  a P-wave  to  the  surface  from  the  depth  of  650  km;  the  corresponding 
feature  for  S can  be  seen  on  Vj  or  W y The  pairs  of  SP  pulses  on  the  trace  correspond  to 
PKIKP  and  pPKIKP;  their  multiple  reflections  are  gradually  attenuated  due  to  finite  Q and 
Internal  reflections  from  discontinuities.  ScS  and  sScS  can  be  easily  seen  on  and  W};  the 
SP  behavior  of  these  two  functions  is  very  similar,  as  could  be  expected  for  Sy  and  S^  at  nearly 
vertical  incidence. 

Once  the  library  of  functions  U,  V,  and  W has  been  prepared,  they  can  be  used  to  compute 
rapidly  synthetic  seismograms  for  an  arbitrary  source  mechanism  and  location  of  the  receiver. 
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The  complete  expressions  are  rather  tedious  and  will  be  omitted  here;  they  are  analogous  to 

a 

Kqs.  (2.1.30)  and  (2.1.31)  of  Gilbert  and  Dziewonski,  and  represent  products  and  sums  of 


functions  IK  f , 
derivatives: 


and  W . . 


with  the  appropriate  spherical  harmonics  of  degree  1 and  their 
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u.(a,e,<M)  = £ Mk(t>  * X Aikf  [<ro.a*t)'(e’0)1  ‘ 

k=l  1=0 


Finally,  after  summation  over  t: 
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ui(r,t)  = ^ Mk(t)  * Bik(r,t* 

k=l 


The  functions  B^(r,t)  represent,  in  effect,  the  Green's  functions  for  each  of  the  six  in- 
dependent elements  of  the  moment  rate  tensor  associated  with  a point  source. 

Figures  11-20  through  11-22  are  record  sections  from  10°  to  180°  computed  for  the  vertical, 
longitudinal,  and  transverse  components  of  an  instrument  with  a WWSSN  response.  The  source 
mechanism  used  was  the  fault -plane  solution  of  Mendiguren^  for  the  Colombian  earthquake  of 
31  July  1970.  A number  of  commonly  observed  body-wave  phases  are  identified  in  the  figures; 
the  effect  of  antipodal  focusing  at  180°  is  quite  spectacular.  Computation  of  54  seismograms 
shown,  each  containing  contributions  of  nearly  6000  normal  modes,  required  only  15  min.  on 
a minicomputer. 

If  the  source  mechanism  can  be  represented  by  a point  source  in  time,  then  the  convolution 

operation  in  the  equation  above  is  replaced  by  multiplication,  and  the  inverse  problem  is  linear 

in  the  time  domain  (it  is  always  linear  in  the  frequency  domain).  Examples  of  application  of  this 
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approach  have  been  demonstrated  by  Stump  and  Johnson  for  sources  imbedded  in  an  elastic 
half-space.  The  following  experiment  is  designed  to  show  the  stability  of  the  inverse  solution 
in  the  presence  of  noise  for  a realistic,  spherically  symmetric  earth  model  with  anelastic 
attenuation. 

Mendiguren's  mechanism  of  the  Colombian  earthquake  (in  equivalent  moment  tensor  rep- 

* Q 

resentation;  cf.  Gilbert  and  Dziewonski,  Eq.  (3.2.1)]  was  used  to  compute  three  3-component 
synthetic  spectra  at  an  epicentral  distance  of  90°  and  azimuths  of  0°,  120°,  and  240°.  The 
inverse  problem  for  the  moment  tensor  was  then  solved  for  various  levels  of  random  noise 
added  to  the  synthetic  spectra. 

Table  II-2  lists  the  results  of  inversion,  with  standard  errors  of  the  elements  of  the  moment 
tensor,  for  different  noise  levels.  A noise  level  of  100  percent  means  that  the  rms  spectral 
amplitude  of  the  noise  is  equal  to  the  rms  amplitude  of  the  signal,  within  the  frequency  band  of 
interest.  The  solution  is  very  close  to  the  exact  values,  even  with  200-percent  noise.  The 
statistically  estimated  standard  errors  are  in  good  agreement  with  the  level  of  differences  be- 
tween computed  and  exact  values.  This  also  applies  to  the  isotropic  component  of  the  moment 
tensor,  whose  value  should  be  zero  for  the  double  couple,  of  course. 

An  explanation  of  the  stability  of  the  solution  is  that  we  have  used  seismograms  that  contain 
body  waves  leaving  the  focus  at  great  many  angles  of  incidence.  Thus,  even  though  we  used  only 
records  from  three  different  sites,  the  sampling  of  the  focal  sphere  was  quite  thorough.  Also, 


It  appears  that  one  could  obtain  a reliable  solution  with  less  information  than  has  been  used  in 
this  particular  example.  This  would  allow  for  exclusion  from  the  analysis  portions  of  the  seis- 
mograms that  are  severely  contaminated  by  the  effect  of  lateral  heterogeneities  — the  funda- 
mental mode  surface  waves,  for  example.  Then,  with  the  source  mechanism  known,  the  phase 
and  group  velocities  could  be  measured  for  the  particular  path.  Such  an  iterative  procedure 
might  represent  a valuable  tool  in  the  studies  on  the  characteristics  of  surface-wave  propagation. 

A.  M.  Dziewonski 

E.  A MOMENT  TENSOR  FORMULATION  OF  ALEKSEYEV'S 

AND  MIKHAYLENKO’S  ELASTIC-WAVE  METHOD 

To  solve  vertically  inhomogeneous  elastic-wave  problems,  A.  S.  Alekseyev  and  B.  G. 
Mikhaylenko1 4 have  proposed  a different  method  which  combines  separation  of  variables  with 
finite  differences.  They  employ  the  Hankel  transform  to  remove  the  radial  dependence  of  the 
solution  as  in  the  reflectivity  method,  but  solve  the  remaining  problem  in  time  and  depth  by 
finite  differences.  The  complete  solution  must  then  be  reconstructed  by  numerically  integrat- 
ing the  inverse  Hankel  transform.  In  a further  development  (unpublished),  they  have  removed 
that  final  numerical  integration  by  a Hankel  series  technique.  Although  it  is  still  too  early  to 
assess  their  method  in  regards  to  its  competitors  (finite  difference,  finite  element,  reflectivity, 
and  Cagnaird-de  Hoop),  it  is  worth  noting  that  it  allows  only  simple  source  models.  Here  we 
attempt  to  reformulate  Alekseyev's  and  Mikhaylenko's  equations  in  a moment  tensor  form  which 
will  allow  any  point  source  mechanism. 

The  equations  of  motion  in  cylindrical  coordinates  can  be  written  as  a Hankel  transform 
+ « 
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which,  in  vector  notation,  reads: 
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where  the  C matrix  and  C vector  are 
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The  closure-completeness  relation  is 
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If  this  form  of  the  solution  is  substituted  into  the  elastic  equation  of  motion  where  the  elastic 
properties  depend  only  on  depth,  the  equations  reduce  to  the  Alekseyev  and  Mikhaylenko  form. 
Strictly  speaking,  however,  the  right-hand-side  force  density  vector  must  be  similarly  trans- 
formed for  a complete  simplification  to  occur. 

In  the,  now  common,  moment  tensor  representation,  the  force  density  vector  is 

written: 


F = — V . [M6(r  - r )]  = 


+ 2 /-•“  — 

£ ^ kdk  Cm(kr,<>)  Hm(k,z,t) 

m=  -2  ° 


where 


Hm(k.z.t)  , 


d(k,z,t,m) 

e(k,z,t,m) 

f(k,z,t,m) 


VM.t)  ^ 


X C*(kr\4>)  V • [M6(r  - r )] 


m 


If  the  moment  tensor  is  expressed  in  north,  east,  and  down  components  as  before,  then  explicit 
expressions  for  the  H vector  are  obtained: 
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Putting  everything  together  gives  the  desired  result: 


which  is  a moment  tensor  form  of  Alekseyev's  and  Mikhaylenko's  result. 

D.  W.  McCowan 
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Fig.  II-4.  Convolution  of  synthetic  wave 
form  with  inverse  filter. 
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Fig.  II-5.  Sample  convolution  of  real  wave  forms  with  inverse  filter 
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Fig.II-6.  Double  couples  superimposed  on  focal  sphere- 
project  of  lower  hemisphere. 
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Fig.II-7.  LP  P-wave  amplitude  from  Cannikin  and  apparent  double-couple 
component:  only  azimuthal  amplitude  variation  is  illustrated,  whereas 
double-couple  component  is  projected  in  equal-area  sense  on  lower  half  of 
focal  sphere. 


Fig.  11-8.  LP  vs  SP  amplitude  variation  for  Aleutian  intermediate -depth 
earthquake,  28  February  1970. 
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Fig.  II- 1 0.  LP  vs  SP  amplitude  variation  for  Borrego  Mountain,  California 
earthquake,  9 April  1968. 
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Fig.II-11.  LP  P-wave  amplitude  from  Novaya  Zemlya  explosion, 
14  October  1970,  and  apparent  double-couple  component. 


Fig.  11-12.  Fault-plane  solution  for  Borrego  Mountain  earthquake.  Equal-area 
projection  of  lower  focal  hemisphere  with  C = compression  and  D = dilation. 
LI  norm  estimate  with  2.7-percent  monopole,  8.3-percent  dipole,  and  89.0-percent 
double  couple. 


Fig  11-13.  Fault-plane  solution  for  Borrego  Mountain  earthquake.  L2  norm 
estimate  with  0.8-percent  dipole  and  99.2-percent  double  couple. 
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Fig.  11-14.  Fault-plane  solution  for  Borrego  Mountain  earthquake.  LI 
norm  estimate  with  6.1-percent  monopole,  16.5-percent  dipole,  and 
77.3-percent  double  couple.  One  data  point  has  been  artificially  ampli- 
fied by  a factor-of-10. 


Fig.  11-15.  Fault-plane  solution  for  Borrego  Mountain  earthquake.  L2 
norm  estimate  with  6.4-percent  dipole  and  93.6-percent  double  couple. 
One  data  point  has  been  artificially  amplified  by  a factor-of-10. 


Fig.  11-16.  Fault-plane  solution  for  Fiji -Tong a earthquake,  9 October  1967. 
LI  norm  estimate  with  17.7-percent  monopole,  23.8-percent  dipole,  and 
58.5-percent  double  couple;  28  LP  amplitudes  and  13  LP  polarities  used. 


Fig.  11-17.  Fault-plane  solution  for  Fiji-Tonga  earthquake.  L2 
norm  estimate  with  4.5-percent  monopole,  7.0-percent  dipole, 
and  87.5-percent  double  couple;  28  LP  amplitudes  used. 
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Fig.  11-18.  Fault-plane  solution  for  Fiji-Tonga  earthquake.  LI 
norm  estimate  with  8.9-percent  monopole,  13.9-percent  dipole, 
and  77. 2 -percent  double  couple;  28  LP  amplitudes,  13  LP  polar- 
ities, and  9 SP  polarities  used. 
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Fig.  11-19.  Excitation  functions  U.  V,  and  W for  source  depth  of  650  km 
and  angular  order  1 = 4.  Elastic  earth  model  1066 B (Ref.  9),  anelastic 
model  FSQMK  (Refs.  5 and  11).  Modes  from  QS4  to  ^gS4were  used  to  con- 
struct functions  U and  V;  functions  W were  obtained  by  superposition  of 
toroidal  modes  from  QT4  to  2iT4'  Note  that  on  U^.  SP  pulses  are  equiv- 
alent to  PKIKP  and  pPKIKP  and  their  multiple  reflections  — ScS  and  sScS 
on  V3  and  W^.  Abrupt  steps  in  Uj,  V3>  and  at  times  of  roughly  1 and 

2 min.  after  origin  time  correspond  to  arrivals  of  P-  and  S-waves,  re- 
spectively, traveling  vertically  from  source  to  surface  of  earth. 
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Fig.  11-20.  Record  sections  of  synthetic  seismograms  for  vertical  component, 
recorded  by  a simulated  WWSSN  instrument,  for  source  depth  of  650  km. 
Double -couple  source  mechanism  determined  by  Mendiguren1^  for  Colombian 
earthquake  of  31  July  1970.  Azimuth  for  all  receivers  is  0°.  Note  amplifica- 
tion of  signal  due  to  convergence  of  waves  at  180°. 


III.  MISCELLANEOUS  STUDIES 


A.  DETECTION  CHARACTERISTICS  OF  THE  SRO  STATIONS 

Several  published  reports1'^  have  described  the  performance  of  the  Seismic  Research 
Observatories  (SROs).  We  present  here  an  updated  evaluation  of  detection  performance,  includ- 
ing that  for  four  new  sites  (ZOBO,  Bolivia;  MAJO,  Japan;  KAAO,  Afghanistan;  and  CHTO, 
Thailand).  For  the  first  86  days  of  1978,  we  have  compared  the  short-period  (SP)  detections 
recorded  at  ten  SRO/ASRO  stations  with  the  SDAC  Bulletin;  both  these  sources  were  accessed 
via  the  Datacomputer. 

For  each  station,  theoretical  P and  PkP  arrival  times  were  computed  for  all  the  events  in 
the  SDAC  Bulletin.  An  event  was  assumed  to  have  been  detected  if  this  arrival  time  lay  within 
the  time  interval  saved.  Since  for  each  SRO/ASRO  detection  the  data  are  saved  starting  10  to 
30  sec  prior  to  the  detection  time,  and  for  several  tens  of  seconds  afterward,  small  errors  in 
origin  and  travel  times  and  location  will  not  affect  the  results.  To  account  for  subsequent  de- 
tection of  later  phases  from  the  same  event,  we  have  also  included  as  "associated"  all  detec- 
tions within  the  subsequent  15  min.  of  a detected  P/PkP  arrival.  Such  secondary-phase  asso- 
ciations accounted  for  2 5 to  50  percent  of  all  associated  detections. 

In  Table  III-l  we  give,  for  each  station  studied,  the  time  period  used  (station  downtimes 
were  not  included),  the  percentage  of  the  total  time  interval  for  which  the  detector  was  on,  the 


TABLE  111-1 

DETECTION  CHARACTERISTICS  OF  TEN  SRO/ASRO  STATIONS 

Period* 

Percent  of  Time 

Number 

Associatedt 

50- Percent 
Detection  Thrediold? 
(mb) 

Station 

Studied 

Detector  On 

of  Detections 

(percent) 

ANMO 

1-86 

3.9 

2058 

50.0 

4.6 

CHTO 

1 -86 

7.0 

4692 

35.3 

4.4 

CTAO 

1 -86 

3.1 

1606 

24.9 

5.2 

GUMO 

1-86 

9.2 

5530 

7.6 

5.8 

KAAO 

1-86 

9.2 

3500 

32.4 

4.4 

MAIO 

60-86 

9.1 

687 

66.0 

4.5 

MAJO 

1-61 

2.9 

630 

41.6 

5.0 

NWAO 

1-85 

3.5 

709 

30.6 

5.7 

TATO 

1-85 

2.4 

1448 

9.3 

5.8 

ZOBO 

1-41 

L 

5.3 

1200 

20.8 



4.8 

* Expressed  at  Julian  days  (1978). 
t See  text  for  <moc lotion  scheme. 

t Estimated  from  data  shown  in  Figs . Ill— 1(a)  through  (j). 
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number  of  detections,  and  the  percentage  of  detections  that  were  associated  according  to  the 
scheme  described  above.  As  has  been  noted  previously/’2  the  detection  performance  of  island 
stations  GUMO  and  TATO  appear  to  be  very  poor. 

At  each  station  we  have  evaluated  the  detection  characteristics  for  events  within  a distance 
of  100°:  the  results  are  shown  in  Figs.  Ill- 1 (a ) through  (j).  In  the  last  column  of  Table  III  - 1 we 
give  the  approximate  50-percent  detection  threshold  inferred  from  the  data  shown  in  these  fig- 
ures. Stations  CHTO,  KAAO,  MAIO,  and  ANMO  have  the  lowest  detection  thresholds,  all  below 
4.6.  Note  that  GUMO,  with  more  total  detections  than  CHTO,  has  an  extremely  high  detec- 

ition  threshold  at  5.8  compared  with  4.4  for  the  latter:  a very  high  proportion  of  the  detections 
for  GUMO  (and  also  NWAO  and  TATO)  are  probably  false  alarms. 

In  an  attempt  to  determine  the  "true"  false-alarm  rate,  we  have  retrieved  from  the  Data-  ! 

computer  the  236  SP  seismograms  corresponding  to  detections  at  ANMO  from  20-31  Decem- 
ber 1976.  Using  the  associated  scheme  described  above,  only  70  of  these  could  be  linked  to 
P/PkP  arrivals,  and  16  to  secondary  phases  of  these  events;  the  50-percent  detection  threshold 
for  this  time  interval  was  ~4.5,  which  is  close  to  that  determined  for  the  longer  time  interval 
[see  Fig.  III-l(a)  and  Table  III-1J. 

Each  of  the  SP  detections  was  bandpass  filtered  over  the  frequency  band  0.5  to  2.0  Hz  to 
simulate  the  data  run  through  the  on-site  detection  processor.2  We  studied  each  of  the  150  un- 
associated detections  to  see  whether  an  analyst  would  have  picked  a phase  in  the  detected  interval.  i 

Four  of  the  "detections"  were,  in  fact,  saved  calibration  pulses.  Of  the  remaining  146  detections, 
inspection  of  the  data  showed  that  an  analyst  as  unexperienced  as  this  author  would  have  picked 
101  as  clear  seismic  phase  arrivals.  For  the  remaining  45,  it  was  difficult  to  see  exactly  what 
triggered  the  detection  algorithm:  this  "false-alarm"  rate  of  19  percent  is  quite  acceptable, 
however,  and  might  even  be  considered  too  low  for  optimum  detector  performance. 

For  the  101  phase  arrivals  selected  above,  we  returned  to  the  raw  (i.e.,  not  bandpass  fil- 
tered) SP  data  and  re-inspected  the  traces.  From  their  high  frequency  character  (dominant 
frequency  >2.5  Hz),  57  of  these  could  be  reasonably  assumed  to  result  from  local  events  at  dis- 
tances of  probably  less  than  10*.  Such  a distance  from  Albuquerque  includes  Southern  California 
and  the  northern  portion  of  the  Gulf  of  California  and  Baja  California.  No  local  Bulletin  for  this 
region  was  available  to  verify  this  hypothesis. 

With  the  exception  of  Australian  stations  NWAO  and  CTAO,  all  the  SRO/ASRO  stations  are 
situated  within,  or  close  to,  regions  of  probably  higher  seismic  activity  than  ANMO.  This  will 
certainly  be  true  for  such  stations  as  CHTO,  GUMO,  KAAO,  MAIO,  MAJO,  and  TATO.  It  is 
thus  likely  that  the  "true"  false-alarm  rates,  i.e.,  percentage  of  detections  which  would  not  be 
picked  by  an  analyst,  are  extremely  low.  A possible  exception  to  this  is  the  notorious  perfor- 
mance of  GUMO,  at  which  the  detector  was  switched  on  for  9.2  percent  of  an  86 -day  interval, 
but  which  detected  only  16  percent  of  events  of  m^  > 4.5  at  distances  of  less  than  100°! 

R.  G.  North 

B.  A MULTIPLE-FREQUENCY  BAND  SHORT-PERIOD  SIGNAL  DETECTOR 

In  Sec.  A above,  we  described  the  results  of  an  evaluation  of  the  detection  performance  of 
the  SRO/ASRO  stations.  It  was  noted  there  that  the  false-alarm  rate  of  the  detector  presently 
operating  at  the  SRO  sites  was  probably  at  least  as  low  as  25  percent,  though  many  of  the  de- 
tections, particularly  those  which  could  be  recognized  as  being  from  close  events  by  their 
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frequency  content,  could  not  be  associated  with  Bulletin  events.  We  present  here  the  results  of 
an  attempt  to  separate  local  and  teleseismic  events  during  the  detection  process. 

Much  research  has  been  carried  out  in  the  past  on  the  design  of  signal  detectors  suitable 
for  seismic  data.  Most  of  these,  however,  have  been  primarily  designed  for  seismic  arrays 
rather  than  the  single  stations  which  are  now  superseding  the  arrays.  All  suggested  detectors 
measure  a particular  property  (such  as  power)  continuously,  and  compare  its  instantaneous 
value  with  a long-term  average  of  this  property,  often  described  simply  by  its  statistical  char- 
acteristics such  as  mean  and  variance.  This  long-term  average  is  either  continuously  or  inter- 
mittently measured,  and  a detection  is  declared  when  the  short-term  instantaneous  value  of  the 
detector  output  exceeds  that  of  the  long  term  by  more  than  a chosen  threshold. 

The  present  SRO  detector  operates  on  a narrowly  (0.5  to  2.0  Hz)  bandpass-filtered  portion 
of  the  SP  filtered  output  of  the  SRO  instrument,  and  measures  power  averaged  over  l-sec  sam- 
ples. Since  only  detections  are  saved,  we  are  generally  unable  to  properly  evaluate  the  per- 
formance of  the  detector  since  we  do  not,  for  example,  know  what  the  character  of  the  "missed" 
signals  (as  predicted  by  Bulletin  location,  for  example)  was  and  precisely  why  the  detector  failed. 
However,  during  the  latter  half  of  1977,  a SRO  instrument  at  ANMO  (Albuquerque)  was  inter- 
mittently run  in  a broadband  mode,  i.e.,  without  separation  into  SP  and  LP  passbands.  We  have 
selected  for  study  a 5f-hr  segment  of  this  continuously  recorded  (20 -Hz)  broadband  data,  start- 
ing at  18:17:51  on  23  December  1977. 

Previous  attempts3  to  utilize  variations  in  spectral  power  as  a detection  method  have  not 
been  shown  to  be  demonstrably  superior  to  the  traditional  technique  of  power  measurement  on 
data  which  have,  as  in  the  case  of  the  SRO  detector,  been  bandpass  filtered  over  frequencies 
assumed  to  be  of  seismic  interest.  Often,  such  spectral  techniques  are  computationally  pro- 
hibitive, such  as  those  which  calculate  spectra  of  moving  or  adjacent  cells3  of  data.  However, 
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a numerically  efficient  scheme  has  been  derived  for  iteratively  computing  spectra,  and  we 
apply  here  a modification  of  this  technique.  The  latter  computed  spectral  power  at  a number  of 
frequencies,  and  used  power  averaged  over  all  frequencies  as  a detection  characteristic.  In 
order  not  only  to  detect  events  but  also  to  decide  whether,  from  their  frequency  content,  they 
are  local  or  teleseismic,  we  have  determined  spectral  power  at  0.5  to  5.0  Hz  in  0.5- Hz  steps  and 
defined  high  (HF)  and  low  (LF)  powers  as 

LF  = | (P(0.5)  + P(1.0)  + P(1.5)  + P(2.0)  J 

HF  = £ [P(2. 5)  + P(3.0)  + P(3. 5)  + P(4.0)] 

where  P(x)  is  th6  spectral  power  at  the  frequency  (x)  in  hertz.  A short-term  average-power 
spectrum  was  determined  from  the  current  60  samples  (3  sec)  of  data,  averaged  over  the  pre- 
vious 60  determinations,  and  the  long-term  power  spectrum  from  the  preceding  600  samples  of 
data,  again  averaged  over  60  samples.  Details  of  the  computation  scheme  may  be  found  in  Ref.  4. 
A detection  was  declared  when  the  ratio  of  either  HF  or  LF  to  its  long-term  average  exceeded  a 
certain  threshold  T.  For  each  detection,  the  maximum  power  in  HF  and  LF  was  measured. 
Figure  III-2  shows  the  amount  of  total  detector  on-time  as  a function  of  threshold  T.  The  de- 
tector has  been  run  on  the  broadband  seismometer  output,  not  on  a filtered  version  of  the  SP  out- 
put channel. 

During  the  time  interval  studied,  the  on-site  detector  at  ANMO  made  10  detections,  totaling 
2500  sec  or  12.9  percent  of  the  total  time  interval.  Of  these,  7 could  be  associated  with 
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SDAC  Bulletin  events  and  3 could  be  readily  seen  to  be  local  events  from  their  high-frequency 
content.  At  T = 3.0,  eleven  individual  detections  were  made  by  the  algorithm  described  above. 

Of  these,  8 corresponded  to  the  ANMO  detections:  the  two  events  missed  were  unfortunately 
teleseismic  Bulletin  events.  At  lower  values  of  the  threshold  T,  the  false-alarm  rate  was  un- 
acceptably high. 

In  Fig.  Ill  - 3 we  show  maximum  values  of  HF  and  LF  for  each  of  the  detections.  Detections 
marked  (A)  were  also  made  by  the  on-site  detector.  Note  that  the  5 teleseismic  events  all  have 
LF/HF  ratios  >1,  and  that  the  3 local  events  have  (LF/HF)  <1.  Although  the  detection  perfor- 
mance of  the  algorithm  is  rather  poor,  it  appears  to  adequately  discriminate  between  teleseismic 
and  local  events.  Similar  results  could,  of  course,  have  been  obtained  by  running  a power  de- 
tector of  the  SRO  type  on  the  outputs  of  two  bandpass  filters  corresponding  to  the  frequency 
ranges  used  in  the  HF  and  LF  determinations.  However,  we  could  have  evaluated  power  spectra 
at  a variety  of  other  frequency  bands  (including,  since  we  are  using  broadband  data,  one  corre- 
sponding to  surface-wave  frequencies)  with  little  more  effort,  and  then  the  iterative  spectral 
power  scheme  would  be  considerably  faster  computationally  than  a series  of  corresponding  band- 
pass filters  on  such  a large  volume  of  data.  Such  a technique  also  has  applications  to  the  compu- 
tation of  spectral  ratios  as  used  in  an  explosion-earthquake  discrimination  contest.  We  propose 
to  further  evaluate  the  above  combined  detection-discrimination  scheme. 

R.  G.  North 


C.  MEASUREMENT  OF  mfe  USING  SRO  DATA 

We  have  commenced  a study  of  the  estimation  of  mb  using  high-quality  SP  waveforms  from 

SRO  and  ASRO  stations.  Several  important  questions  need  to  be  answered.  First,  how  should 

these  waveforms  be  used  in  the  overall  estimation  of  m,  for  an  event  ? Second,  what  are  the 
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characteristics  of  the  SP  waveforms  observed  at  these  sites  from  different  source  areas  ? And 
third,  how  should  we  approach  magnitude  corrections  for  this  network?  We  report  here  some 
preliminary  observations,  based  on  a rather  small  data  set. 

Figures  III— 4 and  1II-5  show  waveforms  observed  at  5 stations  from  two  events  in  the  Aleu- 
tian Islands  (event  details  are  listed  in  the  figure  captions).  Notice  that  Charters  Towers  (CTAO), 
Kabul  (KAAO),  and  Matsushiro  (MAJO)  are  ASRO  sites,  while  the  remainder  are  SROs.  The 
waveforms  are  shown  roughly  equal  in  size,  but  notice  that  there  is  a large  variation  in  the  am- 
plitude scale.  The  approximate  amplitude  of  the  largest  peak-to-peak  value  in  the  first  few 
cycles  is  indicated  in  millimicrons  on  each  waveform.  The  traditional  m^  value,  including  dis- 
tance and  instrument-response  corrections,  is  also  shown  for  each  waveform.  The  mean  SRO 
mb  value  for  each  event  is  close  to  that  quoted  in  the  SDAC  Weekly  Event  Summary  (within  0.1  or 

0.2  m.  units).  Also,  the  standard  deviation  of  the  observations  is  0.3  to  0.4,  which  is  close  to 
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that  usually  observed. 

There  are  clearly  substantial  variations  in  both  amplitude  and  shape  of  waveform.  In  the 
case  of  the  Andreanoff  Islands  event  (Fig.  III-4),  there  is  a substantial  buildup  of  energy  (in  some 
cases  by  an  order-of -magnitude)  over  the  first  20  to  30  sec.  A pP  phase  appears  to  occur  at 
7 to  8 sec  after  the  P-wave  onset,  but  this  seems  generally  larger  than  the  P phase  at  all  azi- 
muths, so  positive  identification  is  not  easy.  The  Fox  Islands  event  (Fig.  Ill- 5 ) is  much  more 
impulsive  at  all  azimuths. 


In  both  cases,  it  is  interesting  to  note  that  large  amplitudes  appear  to  correlate  with  short 
dominant  periods.  In  particular,  Kabul  (KAAO),  Matsushiro  (MAJO),  and  Chiengmai  (CHTO) 
have  dominant  periods  of  close  to  1 sec,  and  larger  amplitudes,  while  the  remaining  stations 
have  dominant  periods  of  up  to  2 sec,  and  much  lower  amplitudes. 

A very  large  difference  is  noted  between  Kabul  (KAAO)  and  Mashhad  (MAIO).  These  sta- 
tions are  relatively  close  together,  and  at  a similar  azimuth  for  events  in  the  Aleutian  Islands, 
yet  the  station  mb  at  Kabul  is  between  0.6  and  1.1  mb  units  higher  than  Mashhad  for  all  events 
studied  in  the  Aleutians  to  date.  The  general  pattern  of  amplitudes  found  for  these  events  is  re- 
peated for  other  events  in  the  Aleutians.  Figure  II1-6  shows  similar  data  for  an  event  in  the 
Kuril  Islands.  Mashhad  is  less  anomalous,  but  large  mb  values  are  still  observed  at  KAAO  and 
CHTO. 

Since  Figs.  III-4  through  III-6  are  typical  of  events  in  the  Aleutians  and  Kurils,  it  would 
appear  that  much  of  the  variation  in  mfa  across  the  network  can  be  attributed  to  systematic  bias, 
rather  than  random  scattering  processes.  North6  was  unable  to  demonstrate  a dependence  of 
receiver  bias  on  source  region.  The  present  results  suggest  a substantial  bias  (perhaps  ap- 
proaching one  magnitude  unit  for  the  Kabul- Mashhad  station  pair)  which  is  consistent  over  a 
source  region  of  the  size  of  the  Aleutian  arc.  Removal  of  this  bias  could  result  in  very  much 
smaller  estimates  of  the  standard  deviation  of  the  random  scattering  component. 

Research  is  continuing  into  the  further  definition  of  these  source  region  dependent  station 
biases. 

J.  C.  Johnston 
M.  A.  Chinnery 

D.  MANTLE  LOVE-WAVE  DISPERSION  FROM  SRO  DATA 

Surface  waves  at  periods  in  excess  of  100  sec  are  an  important  source  of  information  on 
upper-mantle  structure,  since  energy  at  their  wavelengths  (400  to  at  least  2000  km)  is  mainly 
concentrated  in  this  region.  In  particular,  they  are  of  great  use  in  resolving  whether  continent- 
ocean,  and  other,  differences  extend  to  greater  depths  than  shown  by  simpler  crust-uppermost 
mantle  models  deduced  from  shorter-period  surface-wave  dispersion.  However,  such  studies 
have  previously  been  limited  to  recordings  of  the  rare  great  (Mg  > 8.0)  earthquakes.  The  broad 
dynamic  range  of  the  SRO  instrumentation  permits  detection  and  analysis  of  mantle  waves  from 
earthquakes  as  small  as  Mg  = 6.5;  these  are  much  more  frequent  (~15/year  for  1976  to  1978). 
Figure  III-7  shows  an  example  of  mantle  Love  (G)  waves  recorded  at  SNZO  (New  Zealand, 

A = 78*)  from  an  = 7.1  earthquake  occurring  on  26  August  1977  in  the  Southwest  Atlantic 
Ocean.  The  seismogram  f.hown  starts  1000  sec  after  the  event  origin  time  and  is  of  duration 
~7  hr.  The  data  have  been  low-pass  filtered  with  a forward  and  reverse  (phase-free)  Butter- 
worth  filter  with  a 6-dB  point  at  0.01  Hz.  The  second  pass  (G2)  of  the  Love  wave  is  the  largest 
arrival:  subsequent  passes  of  G3  through  G6  are  clearly  visible.  Such  high-quality  data  are 
ideally  suited  to  a study  of  mantle-wave  dispersion  and  Q. 

We  have  analyzed  SRO/ASRO  recordings  of  five  earthquakes  occurring  during  August  to 
December  1977.  The  event  parameters  and  stations  used  are  given  in  Table  III-2  and  the  cir- 
cumferential paths  on  a mercator  projection  are  shown  in  Fig.  III-8.  At  present,  we  have  re- 
stricted ourselves  to  G waves  recorded  on  the  transverse  (after  rotation  into  epicenter  azimuth) 
component;  the  extremely  flat  group- velocity  curve  for  Love  waves  at  mantle -wave  periods  pro- 
duces short-pulse -like  arrivals  as  shown  in  Fig.  III-7.  Each  transverse  component  record  was 
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TABLE  111-2 

EVENT-STATION  PAIRS  USED 


Date 

(1977) 

Origin  Time 

Latitude 

Longitude 

Region 

M 

s 

Stations  Used 

19  August 

06:08:55.2 

11 .1 °S 

1 18.5  “E 

Sumbawa 

7.9 

AN  MO,  SNZO,  CTAO, 
ZOBO,  CHTO 

10  October 

11:53:53.6 

25 .9  °S 

175. 4°W 

Tonga 

7.2 

ANMO,  GUMO,  TATO, 
ZOBO,  NWAO,  CTAO, 
CHTO 

1 October 

17:26:40.4 

27.9°S 

173, 1°E 

Kermadec 

6.7 

ANMO,  GUMO,  MAIO, 
TATO,  NWAO,  CTAO, 
CHTO 

23  November 

09:26:24.7 

31.0°S 

67.8°W 

Argentina 

7.4 

ANMO,  GUMO,  MAIO, 
CTAO,  KAAO,  SNZO 

28  December 

02:45:36.7 

16.7°N 

40.3°E 

Red  Sea 

6.6 

ANMO,  GUMO,  MAIO, 
TATO,  ZOBO,  NWAO, 
CTAO,  KAAO,  CHTO, 
SNZO 

low-pass  filtered  and  reduced  to  a 16-sec  sampling  rate.  Subsequent  passes  of  the  G wave 
(G2  through  G4)  were  windowed  between  group  velocities  of  4.25  to  4.50  km/sec.  We  measured 

7 

phase  velocities  C for  each  record  using  the  classical  formula 


C(T) 


6A 


(5t  + T(|f  + N + i)J 


from  pairs  (G2,  G4)  and  (G3,  G5),  where  6A  is  the  length  of  the  circumferential  path  (corrected 
for  ellipticity),  fit  is  the  difference  in  window  start  times,  and  60  is  the  phase  difference  at  the 
chosen  period  T.  The  factor  of  j accounts  for  the  polar  phase  shift,  and  N is  an  integer  chosen 
to  provide  a smooth  fit  to  free  oscillation  data  at  the  longer  periods.  Such  a differential  mea- 
surement eliminates  both  source  and  instrumental  effects.  The  phase  difference  60  was  deter- 
mined from  the  cross-correlogram  of  the  windowed  seismogram,  as  this  has  been  shown  to 

g 

provide  better  results  than  differences  of  individual  phase  spectra. 

For  each  seismogram,  phase  velocity  was  determined  for  (G4,  G2)  and  (G5,  G3  I pairs,  thus 
providing  a good  check  upor  the  accuracy  of  the  results  over  a single  path.  We  rejected  paths 
for  which  the  phase  velocities  determined  from  (G2,  G4)  and  (G3,  G5)  differed  by  more  than 
0.01  km/sec.  In  Fig.  III-9(a)  we  show  the  results  obtained  for  the  paths  given  in  Table  IH-2  and 
Fig.  Ill -8.  The  results  are  shown  as  deviations  from  the  dispersion  predicted  by  the  oceanic 
model  5.08M  (Ref.  9).  Each  datum  is  the  mean  of  a number  (given)  of  observation.  The  bars 
indicate  ±1  standard  deviation.  Shown  below  this  [Fig.  III-9(b)J,  for  comparison,  are  results 
obtained1®  from  recordings  of  the  great  Alaskan  and  Kuril  earthquakes  of  1964  and  1953.  re- 
spectively, expressed  in  the  same  manner.  We  do  not  attach  any  significance  to  the  deviation 
from  5.08M,  but  wish  rather  to  point  out  the  small  standard  deviation  obtained  using  SRO  data 
from  events  as  small  as  Mfl  = 6.6.  We  propose  to  continue  this  study  (some  45  events  of 
Ms  > 6.5  have  occurred  during  1976-1978)  for  many  more  paths,  including  Rayleigh-wave  data. 
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We  believe  that  enough  results  can  be  obtained  to  permit  a detailed  regionalization  of  phase 
delays  over  several  hundred  paths  in  order  to  contour  phase  velocity  at  a variety  of  periods  and 
(by  implication)  depths  within  the  upper  mantle. 

R.  G.  North 

E.  LITHOSPHERIC  STRUCTURE  OF  THE  WALVIS  RIDGE 
FROM  RAYLEIGH-WAVE  DISPERSION 

{ The  Walvis  ridge  is  a linear  topographic  feature  located  in  the  South  Atlantic  Ocean  (see 

Fig.  Ill  - 10)  which  appears  to  have  no  natural  seismicity.11  The  origin  of  an  aseismic  ridge  like 

the  Walvis  ridge  is  usually  attributed  to  hot-spot  plume  activity  at  mid-ocean  spreading  cen- 
12 

ters.  While  their  kinematic  behavior  agrees  with  plate  tectonic  theory,  the  deep  structure  of 
aseismic  ridges  is  poorly  understood.  Only  one  reversed  marine  refraction  line  exists  sam- 
pling upper-mantle  depths.  The  interpretation  indicates  a 16-km-thick  oceanic  crust  over  a 
low-velocity  (V  = 7.93  km/sec)  mantle.11  The  small  free-air  gravity  anomaly  across  aseismic 
ridges  has  been  used  to  infer  Airy-type  crustal  thickening.1^’19  To  provide  additional  con- 
straints on  the  structure  of  aseismic  ridges,  we  present  here  the  results  from  Rayleigh-wave 
group- velocity  dispersion  analyses.  Unfortunately,  due  to  geographical  limitations,  only  two 
suitable  South  Atlantic  events  were  found  having  nearly  identical  epicenters  (19  April  1968).  So, 
our  results  are  limited  to  just  the  Walvis  ridge. 

The  long-period  vertical  seismograms  from  WWSSN  stations  SDB  and  WIN  were  digitized, 

and  group  velocities  were  calculated  by  the  multiple-filter  technique.10  The  instrument  group 

17  18 

delay  was  removed  by  the  usual  method.  ' Since  the  station  distribution  in  the  South  Atlantic 
is  poor,  no  fault-plane  solutions  could  be  found  and  no  phase  velocities  could  be  measured.  The 
resulting  Rayleigh-wave  group-velocity  dispersion  curves  for  the  on-ridge  (SDB)  and  off-ridge 
(WIN)  crustal  paths  are  compared  with  the  20-  to  50-million-year-old  pure-path  result  of 
Forsyth19  in  Fig.  HI-11. 

Significant  differences  exist  between  the  on-ridge  and  off-ridge  dispersion,  and  between 

both  of  these  and  that  for  a normal  oceanic  lithosphere.  Reasonable  values  of  epicenter  mis- 

iocation,  lithospheric  anisotropy,  and  lateral  refraction  cannot  explain  the  differences.  Station 

timing  was  checked  by  measuring  the  P arrival  time  on  the  SP  record,  and  was  found  not  to 

exceed  1 sec.  This  implies  to  us  that  the  anomalous  dispersion  is  real,  and  that  a deficiency  of 

high-velocity  material  exists  under  the  Walvis  ridge. 

The  dispersion  data  were  inverted  for  the  shear-wave  velocity  structure  using  the  formalism 

of  Wiggins20  as  adapted  by  Rodl  et  al.21  Since  inverting  group-velocity  dispersion  is  a highly 

nonlinear  problem  and  also  because  phase -velocity  information  was  not  available  to  constrain 

the  inversion,  we  searched  for  the  smallest  changes  from  a standard  model  that  fit  the  data. 

19 

The  actual  starting  model  we  used  was  that  of  Forsyth.  Final  models  for  both  South  Atlantic 
paths  are  compared  with  the  20-  to  50-million-year-old  pure-path  result  in  Fig.  Ill- 12.  The  re- 
solving kernels  for  both  models  are  compact,  so  the  parameters  shown  are  required  by  the  data. 
Prominent  features  of  the  on-ridge  model  are  a 12.5-  ± 3-km-thick  crust  and  a low-velocity 
(Vg  = 4.26  km/sec)  upper-mantle  lid  to  depths  of  40  to  50  km.  The  off-ridge  model  has  a normal 
oceanic  crustal  structure  with  a 4.40-km/sec  upper  mantle  to  depths  of  50  km.  No  significant 
differences  from  normal  oceanic  lithosphere  exist  below  50  km  in  either  model. 

Rayleigh  waves  are  sensitive  to  compressional-wave  velocity  and  density  as  well  as  to  shear- 
wave  velocity.  In  order  to  evaluate  the  effect  of  neglecting  these  variables,  the  pertinent  partial 


derivatives  of  group  velocity  were  investigated.  Compressional-wave  velocity  has  no  significant 
effect  below  the  upper  crust  and  therefore  can  be  safely  neglected.  Density  can  affect  the  results 
into  the  upper  mantle  and  must  be  taken  into  account.  Perturbing  the  upper-mantle  density  from 
5.40  to  3.20  gm/cc  produces  a 0.06-km/sec  rise  in  the  upper-mantle  lid  velocity  and  a very  slight 
reduction  below  this  depth.  In  either  case,  the  resulting  shear  velocities  shown  in  Fig.  III-12  are 
significantly  lower  than  those  of  normal  oceanic  lithosphere. 

I’he  velocities  presented  here  are  averages  for  the  entire  paths  shown  in  Fig.  Ill- 1 0.  The 
differences  between  the  on-ridge  and  off-ridge  results  may  be  due  to  the  Rayleigh  waves  travers- 
ing a path  which  is  half  over  the  ridge  and  half  over  the,  presumably  normal.  Cape  Basin.  The 
available  data  are  not  capable  of  distinguishing  this  effect. 

Possible  explanations  for  the  low-velocity  mantle  in  our  model  include  partial  melting  and 
compositional  changes.  The  former  can  be  eliminated  on  the  basis  of  the  measured  heat  flow, 
which  is  only  slightly  higher  than  that  in  the  surrounding  ocean  basins.22  The  second  possibility 
is  attractive  due  to  the  unique  isotopic  composition  of  Walvis  ridge  volcanics23'24  which  could 
indicate  an  anomalous  source  region.  Iron  enrichment  in  the  mantle  can  easily  produce  the  re- 
quired reduction  in  shear  velocity  with  no  large  density  change.25  Further  data,  preferably 
compressional-wave  velocity  measurements,  are  required  to  settle  the  issue. 

A.  D.  Chavet 
J.  D.  Phillips 
D.  W.  McCowan 

F.  PHASE  VELOCITY  OF  Lg  IN  NORTH  AMERICA 

Maximum-likelihood-array  processing26  has  been  used  to  determine  the  phase  velocity  of 
the  vertical  component  of  Lg  for  several  North  American  continental  paths.  Table  III-3  lists 


TABLE  (11-3 


ARRAY  ANALYSIS  EVENT  PARAMETERS 


Event 

Origin  Time 

Array 

A 

(deg) 

V (Group) 
(krrv/sec ) 

V (Phare) 
(kn/sec) 

AV  (Phase) 
(kny/tec) 

Frequency 

(Hz) 

GAS 

10  December  1967 
19:30:00.1 

LASA 

10.0 

3.5 

1 

4.4 

0.18 

1.1 

MINN 

9 July  1975 
14:54:16.3 

YKC  ; 

1 

20.0 

3.5 

4.5 

i 

i 

0.05 

1.0 

N.C. 

30  November  1973 
07:48:38.7 

LASA 

19.7 

3.5 

1 

3.9 

0.08 

1.1 

15  July  1971 
00:24:03.1 

LASA 

! 19.2 

3.5 

4.6 

0.11 

0.9 

the  events  used  in  this  study  and  the  corresponding  measured  phase  velocity  and  precision 
obtained  at  the  listed  array.  The  LASA  B2  subarray  was  used  for  the  Gasbuggy  (GAS)  explo- 
sion and  the  North  Carolina  (N.C.)  and  Queen  Charlotte  Island  (Q.C.)  earthquakes.  The  Yellow- 
knife (YKO)  array  was  used  to  measure  the  phase  velocity  of  the  Minnesota  (MINN)  earthquake. 

In  each  case,  the  measurement  was  made  at  the  time  corresponding  to  an  Lg  group  velocity 
of  3.5  km/sec,  and  at  a frequency  near  1 Hz.  For  every  event  except  GAS,  the  L^  arrival 
was  well  above  the  coda  level  with  its  maximum  amplitude  near  the  above  group  velocity  of 
3.5  km/Bec.  At  the  time  corresponding  to  that  group  velocity  for  the  GAS  data,  there  was  no 
discernible  change  in  the  coda  level;  however,  the  array  processing  clearly  indicated  the  L^ 
arrival.  The  precision  of  the  measurements  was  determined  by  converting  the  3-dI3  half-width 
of  the  corresponding  wavenumber  power  spectrum  peak  to  velocity.  The  results  are  shown  in 
Fig.  III-13.  In  the  previous  SATS,27  the  same  measurement  was  made  for  a Western  Russia 
event  at  a subarray  of  NORSAR.  For  those  data  at  the  same  3.5-km/sec  group-velocity  point, 

the  phase  velocity  was  4.1  km/sec. 

27-29 

As  has  been  noted  by  others,  these  measurements  are  consistent  with  the  behavior  of 
higher-mode  Rayleigh  waves  whose  main  lobe  of  elastic  energy  is  confined  to  the  middle  portion 
of  the  crust.  However,  the  phase  velocities  at  1 Hz  are  significantly  higher  than  those  that 
would  be  obtained  for  a J-B  model  of  the  crust  and  upper  mantle  unless  modified  by  a low-shear- 
velocity  zone  (LVZ)  in  the  upper  mantle  overlaid  by  a high  (4.8-km/sec  or  more)  shear-velocity 
lid.1 2 3 4 5 6^  These  results  support  the  existence  of  a LVZ  with  a high-velocity  lid  in  the  most  stable 
of  continental  regions  (the  Baltic  and  Canadian  Shields),  perhaps  implying  that  a well-defined 
lithosphere-asthenosphere  boundary  persists  continuously  around  the  globe.  Our  data  also 
tend  to  substantiate  that  these  waves  propagate  coherently  along  apparently  high-Q  paths  over 
tectonic  environments  ranging  from  Shield  to  Basin  and  Range  regimes.  On  this  basis,  coher- 
ent processing  of  the  phase  should  aid  discrimination  and  detection  capability  for  continental 
source-receiver  paths.  T.  E. 

T.  J.  Fitch 
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Fig.  Ill-l(a-j).  Detection  performance  of  ten  SRO/ASRO  stations,  mb-frequency  statistics 
for  events  within  100*  (source:  SDAC  Bulletin).  Black  portion  indicates  detection,  white 
indicates  nondetection.  Time  interval  studied  (Julian  date,  1978)  is  given  for  each  station. 
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Fig.  ni-2.  Detector  on-time  as  a function 
of  detection  threshold  T.  Lowest  false- 
alarm  rate  consistent  with  best  perfor- 
mance obtained  with  T = 3.0. 
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Fig.  ID-3.  High  (HF)  and  low  (LF)  frequency 
spectral  power  averages  measured  at  peak 
power  (larger  of  HF  and  LF)  for  each  of  11 
detections  obtained.  Detections  also  made 
by  SRO  detector  are  designated  (A)  - mag- 
nitude and  location  of  associated  detections 
are  also  given.  Note  that  teleseismic  ar- 
rivals lie  above  line  of  HF  = LF. 


Fig.  III-7.  Low-pass  filtered  transverse  component  recording  at  SNZO 
of  earthquake  in  southwest  Atlantic  Ocean.  Seismogram  starts  4000  sec 
after  origin  time,  is  ~ 7 hr  in  duration.  Love  waves  G2  through  G6  shown. 
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Fig.  DI-8.  Circumferential  paths  for  event-station  pairs  given  in  Table  M-2 
(mercator  projection).  Sampling  is  more  dense  for  Southern  Asia,  Northern 
Africa,  South  America,  and  the  Southern  Pacific.  More-uniform  sampling 
will  be  obtained  by  an  extension  of  this  study. 
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Fig.  Ill— 9.  Love-wave  phase  velocities,  expressed  as  means 
and  standard  deviations,  with  respect  to  model  5.08M  (Ref.  9). 
Present  results  are  shown  in  (a),  those  of  a previous  study!0 
in  (b).  Number  of  observations  at  each  datum  is  given. 
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Fie.  III-10.  Index  bathymetry  chart  of  South  Atlantic  Ocean  showing  Walvts 
ridee  extending  southwestward  from  African  margin  to  mid -Atlantic  ridge. 
Depths  contours  are  in  thousands  of  meters.  Depths  shallower  than  300 m 
are  stippled.  Curved  arcs  show  ray  paths  between  mid-Atlantic  ridge  earth- 
quake epicenters  (cross)  and  WWSSN  stations  SDB  and  WIN. 


IV.  COMPUTER  SYSTEMS  AND  DATA  COLLECTION 


A.  A UNIX  DATA-ANALYSIS  AND  DISPLAY  SYSTEM 

The  imminent  demise  of  the  PDP-7  computer  systems  necessitates  the  transfer  of  a number 
of  software  packages  from  the  PDP-7  system  to  UNIX.  One  of  the  most  notable  of  these  is  the 
Data  Analysis  and  Display  System  (DADS).  The  following  section  is  a brief  description  of  a ver- 
sion of  the  DADS  that  we  currently  are  designing  for  UNIX. 

1.  General  Description 

The  UNIX  DADS  is  a data-management  system  designed  specifically  for  seismic  databases. 
The  system  is  general  enough  to  store  the  large  variety  of  data  and  data  types  needed  in  seismic 
databases.  It  is  also  specific  enough  to  take  advantage  of  some  of  the  efficiencies  made  possible 
by  the  basic  properties  of  seismic  databases.  DADS  makes  no  claim  at  being  a general  data- 
management  system. 

Conceptually,  DADS  may  be  thought  of  as  a two-dimensional  data  matrix.  Matrix  rows  rep- 
resent data  objects  whose  properties  are  stored  in  the  database.  Each  matrix  set  of  values  is 
used  to  delineate  the  properties  of  the  data  objects.  Each  matrix  column  is  defined  as  having 
an  explicit  data  type  such  as  time,  depth,  or  magnitude. 

A typical  DADS  data  matrix  is  a seismic  event  list.  Each  row  in  the  event  list  represents 
some  seismic  event.  The  columns  represent  the  properties  of  the  events  such  as  depth,  time, 
and  magnitude.  Columns  may  also  be  used  to  represent  larger  entities  such  as  lists  of  arrival 
times  for  each  event,  or  even  a list  of  waveforms  for  that  event. 

There  are  a variety  of  commands  which  may  be  used  to  manipulate  DADS  databases.  A 
command  line  interpreter  accepts  a set  of  DADS  commands  and  executes  them.  These  commands 
specify  subsets  or  supersets  of  databases  and  what  is  to  be  done  with  them.  The  user  may  wish 
to  search  the  database  for  a particular  subset  of  data  and  store  it  in  a new  database. 

Other  operations  of  interest  are  the  "table"  command,  which  allows  users  to  display  the 
specified  data  set  in  table  form,  a graph  command,  which  allows  the  user  to  plot  two  variables, 
and  a set  of  commands  to  perform  arithmetic  calculations. 

2.  Implementation:  File  Structure 

The  DADS  implementation  relies  on  a number  of  fundamental  assumptions.  The  most- 
important  assumption  is  that  most  of  the  computation  done  by  the  DADS  software  will  be  in  the 
form  of  database  searches.  The  amount  of  processing  needed  to  generate  a new  textual  display 
of  a database  is  negligible  compared  with  the  processing  required  for  database  searches.  An- 
other assumption  is  the  type  of  terminals  which  will  be  used  as  DADS  consoles.  DADS  is  de- 
signed to  work  well  with  Tektronix  4014  and  VT52  terminals.  It  does  not  work  with  Decwriters 
or  any  other  hardcopy  terminals. 

The  file  structure  used  to  implement  DADS  reflects  the  efficiency  considerations  mentioned 
above.  Each  database  is  a file  directory  which  contains  a collection  of  files  of  three  types.  The 
first  file  type  is  the  header  file  which  contains  a description  of  the  database  or  database  subset, 
along  with  bookkeeping  information.  The  second  file  type  is  the  data  file  which  is  used  to  store 
database  data.  The  last  file  type  is  the  selector  file  that  contains  the  information  used  to  select 
which  portions  of  the  data  belong  to  which  database  subsets. 
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The  data  files  are  organized  by  data  matrix  column.  Each  data  file  is  a one-dimensional 
array  whose  elements  delineate  the  data  objects  in  the  database.  The  constituent  data  elements 
of  a data  object  are  the  elements  with  the  same  index  within  each  data  file  specified  by  the  header 
file.  This  configuration  of  data  files  has  a number  of  advantages  for  processing  the  data.  A 
search  or  scan  of  the  values  of  one  data  descriptor  is  reduced  to  a sequential  scan  of  one  data 
file.  The  addition  or  deletion  of  a data  column  involves  the  simple  addition  or  deletion  of  a file 
to  the  database.  New  entries  in  the  database  are  implemented  by  appending  a new  data  element 
onto  each  data  file. 

The  header  files  along  with  the  data  selector  file  are  the  key  parts  of  the  database  used  to 
organize  the  data.  The  header  file  includes  an  identifying  text  string  of  up  to  512  characters, 
a unique  selector  identifier  (ID),  and  a set  of  data-type  information  for  the  data  column  files. 
There  is  one  master  header  file  which  is  used  as  an  index  of  the  entire  database.  The  other 
header  files  are  used  to  select  database  subsets.  The  master  header  file  always  has  the  same 
name  as  the  database  directory.  The  DADS  is  currently  limited  to  1 database  master  header 
file  and  30  database  subset  header  files. 

The  selector  ID  is  used  in  conjunction  with  the  selector  file  to  identify  which  of  the  data 
objects  belongs  to  a particular  database  subset.  The  selector  file  is  a one -dimensional  array 
similar  to  the  data  files.  Each  element  in  the  array  contains  a key  that  identifies  the  database 
subsets  to  which  the  corresponding  elements  in  the  data  files  belong.  The  corresponding  ele- 
ments in  the  data  files  are  those  elements  with  the  same  file  index  as  the  selector  file  element. 
The  selector  ID  itself  is  a long  integer  with  one  bit  set  to  1.  Each  database  subset  is  assigned 
a different  bit.  A selector  match  occurs  if  the  logical  AND  of  the  selector  ID  and  the  selector 
file  element  is  nonzero.  A scan  of  any  database  subset  will  therefore  begin  with  a scan  of  the 
selector  file  to  see  if  the  data  in  any  particular  row  are  elements  of  that  subset. 

The  data  column  information  in  a header  file  consists  of  a character  string  of  up  to  8 char- 
acters and  an  integer  data  type  code.  The  8 characters  correspond  to  data  file  names  in  the 
database  directory.  The  data  type  code  corresponds  to  the  legal  data  type  codes  such  as  integer, 
real  float,  etc.  which  may  be  found  in  a data  file. 

A sample  database  is  shown  in  Fig.  IV-1.  The  database  is  an  event  list  with  7 data  columns 
of  varying  data  types.  Also  included  are  two  database  subset  groupings  of  three  data  columns 
each.  The  selector  file  indicates  which  rows  in  the  database  correspond  to  which  subsets  of  the 
database.  In  the  example,  00  Indicates  that  a row  is  only  a member  of  the  whole  database  and 
not  in  any  subset.  A 01  indicates  that  a row  is  a member  of  subset  1.  A 10  indicates  that  a 
row  is  a member  of  subset  2.  All  indicates  that  a row  is  a member  of  both  subsets.  The  last 
column  is  a special  data  type.  It  contains  the  names  of  waveform  files  as  elements. 

3.  Implementation:  Command-Line  Interpreter 

A simple  interpretive  language  is  used  to  access  and  manipulate  DADS  databases.  The 
language  has  a limited  syntax.  It  does  not  allow  conditional,  branch,  or  iteration  statements. 
These  features  will  be  added  at  some  future  date.  The  language  does  allow  the  selection  and 
manipulation  of  database  subsets.  It  also  interprets  arithmetic  expressions  written  in  standard 
infix  notation. 

While  any  database  in  the  system  may  be  accessed  with  a command  string,  it  is  more  ef- 
ficient to  specify  the  current  database  of  interest.  This  may  be  done  with  the  command  line: 
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current  events; 

where  "events"  is  the  name  of  a database  header  file  in  the  local  database  directory.  The  cur- 
rent database  can  be  changed  at  any  time.  References  which  do  not  use  the  current  database 
must  include  the  relative  pathname  of  that  database.  These  references  will  be  considerably 
less  efficient.  It  is  also  recommended  that  the  user  move  to  the  directory  in  which  a database 
is  located  before  invoking  the  DADS  on  that  directory.  The  DADS  is  invoked  with  the  shell 
command: 

dads  "databasename"; 

The  manipulation  of  a database  is  accomplished  through  a combination  of  Boolean  expres- 
sions and  assignment  statements.  The  Boolean  expressions  are  used  to  select  a subset  of  the 
database.  The  subset  may  include  any  number  or  combination  of  data  matrix  columns  or  rows. 
The  method  used  for  specifying  rows  is  a comparison  test.  An  example  of  this  is  the  selection 
of  a data  subset  from  an  event  list  database.  If  only  those  events  between  I960  and  1970  with 
magnitude  greater  than  7 are  to  be  selected,  one  would  use  the  following  Boolean  expression: 

magnitude  > 7 fcfc  I960  < year  < 1970; 

If  only  the  depth,  magnitude,  and  year  data  columns  are  required,  one  would  expand  the  state- 
ment to: 

magnitude  > 7 &&  1960  < year  < 1970:depth,  magnitude,  year; 

The  target  of  this  database  subset  could  be  either  a stored  database  subset  or  a display  of 
the  data.  A stored  subset  may  be  created  by  declaring  a subset  database  header  variable  and 
assigning  it  the  result  of  the  Boolean  expression: 

subset  foo; 

foo  a magnitude  > 7 I960  < year  < 1970:depth,  magnitude,  year; 

The  database  or  database  subset  is  displayed  with  the  "table"  command.  This  command 
acts  as  a function  call  and  can  be  used  with  either  a database  name  or  a Boolean  expression  as 
a parameter.  Two  examples  of  this  are: 

table(foo); 

table  (magnitude  > 7 1960  < year  1970:depth,  magnitude,  year); 

In  the  first  example,  the  database  with  name  foo  is  displayed.  In  the  second  example,  the  sub- 
set of  the  current  database  as  specified  by  the  Boolean  expression  is  displayed.  "Table"  can 
also  work  on  individual  data  columns  within  a database.  The  command: 

table(:year); 

is  also  legal  if  year  is  a database  column.  In  this  case,  one  entire  database  column  (or  as  much 
as  will  fit  on  a screen)  is  displayed. 

Individual  database  columns  can  also  be  manipulated  through  interpretive  language.  For 
example,  the  command  lines: 

column  logdepth 
logdepth  = log(depth); 

define  a new  data  column  logdepth  in  the  current  database  and  store  the  log  of  all  the  correspond- 
ing entries  of  the  depth  column.  The  data  type  of  logdepth  will  be  coerced  to  the  same  data  type 


as  the  depth  column.  It  is  also  possible  to  specify  columns  from  a database  other  than  the  cur- 
rent one.  This  is  done  by  specifying  the  relative  pathname  to  each  data  column  involved.  For 
example,  the  command  line: 

. ./events/logdepthll]  = log(depth(  1 )); 

assigns  the  log  of  the  first  entry  of  the  "depth"  column  in  the  current  database  to  the  third'entry 
in  the  "logdepth"  column  which  is  located  in  the  "events"  database.  The  events  database  is  a 
sibling  of  the  current  directory.  This  command  line  also  demonstrates  how  individual  elements 
within  a data  column  can  be  accessed. 

The  arithmetic  expressions  allowed  by  the  interpretive  language  include  the  standard  op- 
erators such  as  + - * / as  well  as  Boolean  operators  and  certain  function  operators.  The  user 
may  declare  variables  in  any  legal  database  datatype  and  use  these  variables  in  the  arithmetic 
expressions.  Data  columns  and  database  header  file  names  are  permissible  data  types  for 
certain  types  of  arithmetic  expressions. 

J • oSX 

B.  UNIX  SIGNAL  DISPLAY  PACKAGE 

During  the  past  six  months,  we  have  made  three  significant  improvements  in  the  UNIX  Sig- 
nal Display  Package.^ 

First,  a cursor  feature  has  been  implemented.  The  user  can  now  call  up  the  Tektronix 
cursor  and  use  it  to  set  named  markers  which  visually  flag  particular  data  points  of  a seismo- 
gram. The  locations  of  these  markers  become  a permanent  part  of  the  database  and  may  be 
used  as  input  to  related  seismic -processing  programs. 

Second,  the  Display  program  has  been  modified  to  utilize  the  Versatek  capability  of  our 
graphics  software.  There  are  two  options:  the  user  may  choose  to  get  a Versatek  copy  of  the 
current  picture  on  his  screen,  or  he  may  choose  to  extend  the  waveforms  off  to  the  right  for  a 
maximum  output  of  9 Versatek  pages  (76  in.). 

Third,  the  speed  of  the  Display  program  has  been  increased  substantially  by  eliminating 
options  that  were  no  longer  needed  and  improving  the  efficiency  of  certain  I/O  operations. 

L.  J.  Turek 

C.  TEKTRONIX  GRAPHICS  SPEEDUP 

The  UNIX  graphics  system  is  currently  too  slow  to  be  used  as  a truly  interactive  system. 
Although  there  are  many  contributing  factors  to  this  problem,  one  of  the  major  problem  areas 
is  the  communication  link  between  the  PDP-11  processor  and  the  Tektronix  graphics  terminal. 
This  link  currently  consists  of  a standard  9600 -baud  teletype  interface.  Analysis  of  the  system 
software  and  hardware  indicates  that  this  communication  link  is  too  slow  to  effectively  service 
the  Tektronix  4014  terminal. 

The  system  software  is  particularly  inefficient  in  the  way  it  handles  character  device  inter- 
faces. While  the  interface  is  acceptable  for  normal  teletype  usage,  the  software  overhead  of 
handling  character  buffers  is  particularly  painful.  The  character  device  Interface  was  not  de- 
signed to  handle  the  large  volume  of  data  required  to  transmit  a graphics  picture  to  the  Tektronix. 


t Seismic  Discrimination  Semiannual  Technical  Summary,  Lincoln  Laboratory,  M.I.T. 
(31  March  1978),  DDC  AD-A057279. 


69 


■NBWMU1  iMi»  iunii  ii 


The  hardware  also  contributes  to  the  problem.  The  Tektronix  4014  is  capable  of  accepting 
graphics  characters  at  rates  approaching  50K  baud  rather  than  the  9600  baud  provided  by  the 
teletype  controller.  The  increased  data  rate  can  be  achieved  by  synchronizing  the  data  flow 
from  the  CPU  with  the  rate  at  which  the  Tektronix  screen  is  painted.  The  synchronization  is 
necessary  because  the  rate  at  which  a Tektronix  can  accept  data  varies  with  the  rate  at  which 
it  draws  vectors.  The  vecto'r  drawing  rate,  in  turn,  varies  with  the  size  of  the  vectors. 

A properly  designed  interface  between  the  PDP-11  and  the  Tektronix  terminal  should  go  a 
long  way  toward  speeding  up  graphics.  The  design  must  include  both  software  and  hardware 
modifications  to  the  PDP-11  to  Tektronix  interf  The  software  character  handling  overhead 
must  be  eliminated,  and  hardware  data  synchr  on  mechanism  must  be  developed. 


1.  Interface  Hardware 

The  hardware  configuration  of  the  graphics  interface  consists  of  a DMA  I/O  device  at  the 
CPU  end,  and  a parallel  interface  with  synchronization  control  at  the  Tektronix  terminal.  The 
DMA  device  has  a software  interface  similar  to  the  DEC  DRlt-b.  Its  command  and  status  reg- 
isters will  perform  basically  the  same  functions  as  the  DRll-b.  In  addition,  the  DMA  device 
has  a FIFO  buffer  of  about  512  bytes.  The  FIFO  is  used  to  store  data  waiting  to  be  transferred 
to  the  Tektronix. 

The  Tektronix  interface  is  the  general-purpose  parallel  interface  021-0109-00.  It  has  the 
necessary  control  logic  to  handle  up  to  8 separate  terminals.  Each  terminal  has  a separate 
address  which  the  CPU  specifies  when  transmitting  data  to  a terminal. 

2.  Interface  Software 

The  interface  software  treats  the  DMA  device  as  a block  I/O  device  rather  than  a character 
I/O  device.  This  allows  for  the  efficient  transfer  of  up  to  512  bytes  directly  from  memory. 
While  these  block  transfers  are  somewhat  wasteful  for  normal  teletype  interactions  with  the 
Tektronix  terminal,  they  are  absolutely  essential  for  efficient  data  transfer  of  graphics  data. 

The  user  program  interface  to  the  Tektronix  terminals  is  totally  unaffected  by  the  software 
and  hardware  modifications. 


3.  Data  Transfers 

A data  transfer  from  the  CPU  to  the  terminal  begins  with  a DMA  transfer  of  a block  of  data 
from  memory  to  the  FIFO.  This  causes  a data-ready  flag  to  be  raised  on  the  Tektronix  inter- 
face. The  first  byte  of  the  data  specifies  which  terminal  is  to  receive  the  data.  The  specified 
terminal  gains  control  of  the  Tektronix  data  interface.  The  terminal  initiates  a data  transfer 
by  strobing  the  FIFO  register  for  each  byte  of  data  to  be  transmitted.  The  strobe  rate  is  gov- 
erned by  the  rate  at  which  the  data  can  be  executed  as  graphics  commands  by  the  terminal.  The 
DMA  resets  the  data-ready  line  after  the  last  byte  in  the  FIFO  has  been  transmitted  to  the  Tek- 
tronix. The  DMA  interface  also  generates  an  interrupt  to  the  CPU  to  indicate  that  the  data  trans 
fer  is  complete.  The  CPU  may  now  send  the  next  data  block  to  the  FIFO  for  transmission  to  a 
Tektronix  terminal. 

A timeout  mechanism  is  included  in  the  DMA  device  for  measuring  the  response  time  of  a 
Tektronix  terminal.  If  the  terminal  does  not  respond  within  1 sec  to  a data-ready  line  or  if  the 
time  between  data  strobe  pulses  exceeds  1 sec,  the  terminal  is  declared  down,  the  FIFO  will  be 


flushed,  and  an  error  interrupt  is  sent  to  the  CPU  to  indicate  a failure  in  the  data  transmission. 
This  mechanism  will  prevent  a dead  terminal  from  tying  up  the  graphics  interface. 

Data  transmission  from  a Tektronix  terminal  to  the  CPU  is  initiated  by  the  terminal,  which 
sends  an  interrupt  on  one  of  8 lines  to  the  DMA  device.  Upon  receiving  an  interrupt,  the  DMA 
device  strobes  the  appropriate  terminal  to  transmit  the  data  to  PDP-11  memory  via  the  DMA 
input  registers. 

The  overall  interface  design  is  relatively  simple  and  could  be  completed  with  only  a modest 
investment  in  hardware  and  software  design  and  implementation.  The  only  significant  difficulty 
will  be  in  building  the  DMA  device  for  the  PDP-11.  Standard  DMA  interfaces  built  by  DEC  or 
other  manufacturers  do  not  include  all  the  control  logic  necessary  for  this  particular  interface. 
One  possible  starting  point  is  the  MDBll-b  general-purpose  DMA  controller  which  provides  a 
generalized  DMA  interface  for  the  PDP-11  UNIBUS,  including  the  necessary  control  and  status 
registers.  This  device  can  be  used  in  conjunction  with  the  MDBU-ww  wirewrap  expansion 
module  which  is  needed  to  hold  the  additional  circuitry  for  the  FIFO  and  the  control  logic  to  be 
used  in  conjunction  with  the  Tektronix  interface.  The  entire  package  makes  possible  a signifi- 
cant speedup  of  the  Tektronix  graphics.  . _ 


D.  DATACOMPUTER  SOFTWARE 

Access  to  the  SRO  raw  data  files  has  been  improved.  In  the  31  March  1978  SATS,  we  de- 
scribed a new  program  (called  DCRISP)  which  outputs  properly  formatted  Datacomputer  short- 
period  (SP)  request  files  by  computing  the  theoretical  arrivals  at  SRO  stations  for  an  input  event 
list.  We  have  rounded  out  this  project  by  adding  a second  program  (called  DCRILP)  which  per- 
forms the  same  function  for  long-period  (LP)  arrivals. 

We  have  also  added  a new  program,  called  WFDAY,  which  will  access  the  Preliminary 
Summary  Waveform  Files  that  are  now  being  created.  These  are  files  that  will  eventually  con- 
tain all  the  waveform  segments  associated  with  a given  event.  SRO  data  are  not  yet  being  stored 
in  these  files,  but  data  from  certain  Alaskan  stations  are  stored  in  these  files  and  nowhere  else. 

WFDAY  has  two  modes  of  operation:  it  will  either  retrieve  all  the  waveform  segments  for 
a given  day,  or  it  will  retrieve  all  the  waveform  segments  stored  for  a specified  event.  The 
resulting  data  will  be  reformatted  and  stored  in  a UNIX  waveform  database,  ready  for  viewing 
with  the  Display  program,  etc.  , , ~ k 


E.  TAPE  DATA  COLLECTION 

Until  recently,  the  Applied  Seismology  Group  has  housed  a library  of  nearly  i?,000  digital 
tape  reels  of  data  from  LASA,  NORSAR,  and  ALPA.  It  has  not  been  possible  to  carry  out  ade- 
quate tape  maintenance  on  a collection  of  this  size.  As  a result,  the  earlier  portions  of  data 
have  been  subject  to  deterioration.  Specific  problems  noted  have  included  surface  deformations 
and  shedding  of  the  oxide  coating,  which  result  in  extensive  data  dropouts,  and  some  complete 
voids. 

As  a result  of  these  problems,  and  the  fact  that  use  of  these  data  has  decreased  substan- 
tially in  the  last  few  years,  we  have  decided  to  make  a significant  reduction  in  the  size  of  the 
data  collection.  With  the  exception  of  data  for  underground  explosions,  all  LASA  data  for  the 
period  1969  through  1969  are  being  removed  from  the  library,  as  are  all  reels  containing  noise 
and  other  special-purpose  Information.  Approximately  5000  of  these  tape  reels  contain  useful 


information.  We  hope  to  identify  an  archiving  agency  which  will  assume  the  storage  and  main- 
tenance of  these  tapes. 

The  collection  vvhich  we  will  maintain  consists  of  LASA,  NORSAR,  and  ALPA  data  for  the 
years  1970  through  1973,  plus  the  explosion -data  set.  This  collection  will  number  just  over 
3000  tape  reels.  The  maintenance  program  will  consist  of  respooling  and  cleaning  the  tapes 
every  year,  and  transferral  to  new  tapes  every  10  years. 

The  array  data  collection  is  now  being  supplemented  by  the  network  day  tapes  from  the  SRO 
network,  which  are  being  received  on  a regular  basis. 
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